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ABSTRACT 

A  quasi-one  dimensional,  nonequilibrium,  streamtube  model  and  finite 
difference  computer  code  has  been  developed  to  study  nonisentropic 
supercritical,  real  gas  flows.  A  system  of  unsteady  Euler  equations  is 
integrated  to  a  steady  state  solution  utilizing  MacCormack’s 
explicit/implicit  numerical  scheme.  Implicit  characteristic  boundary 
conditions  are  employed  at  the  subsonic  inlet  to  enable  automatic  solution 
of  Inlet  conditions  consistent  with  the  critical  mass  flow  rate.  A 
numerical  damping  technique  for  the  explicity  evaluated  reaction  rate  has 
been  developed  such  that  nonequllibrluia  flows  with  large  reaction  rate 
coefficients  have  been  solved  numerically  with  significant  reduction  in 
computation  time. 
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NOMENCLATURE 

streamtube  cross-section  area,  also  symbol  for  atomic  species 
speed  of  sound 

frozen  specific  heat  at  constant  volume 
specific  internal  energy 
specific  enthalpy,  also  Flank's  constant 
cut  off  coefficient 
equilibrium  constant 

forward  rate  constant  for  reaction  A2  +  A^2A  +  A 

forward  rate  constant  for  reaction  A2  +  A2^  2A  +  A2 

Boltzmann's  constant 

characteristic  length 

molecular  weight  for  atomic  species 

mass  flow  rate 

pressure 

partition  function 

specific  gas  constant  of  the  mixture 

entropy 

temperature 

time 

velocity  in  the  streamwise  direction 
volume 

atomic  species  equation  source  term 
streamwise  coordinate 
cross-section  coordinate 


a 


mass  fraction  of  atomic  species 
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6  proportionality  constant  for  integration  step  size  baaed  on 

chemical  relaxation  time 

y  ratio  of  specific  heats 

6  characteristic  temperature 

p  density 

T  charac t  eris  t ic  t  ime 

^  reaction  rate  coefficient 


Subscripts 

c  critical 

d  dissociation 

el  electronic 

F  frozen 

i  z-mesh  point  index 

rot  rotational 

t  throat 

tr  translational 

vib  vibrat lonal 

a  differentiation  with  respect  to  a 

p  differentiation  with  respect  to  p 

0  stagnation 

Superscripts 

n  time  level  index 

nf 

* 


temperature  exponent 
equilibrium 
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Chapter  I 

INTRODUCTION  AND  LITERATURE  SURVEY 


In  the  past  few  decades  Interest  In  high  temperature  aerodynamics  has 
Increased  dramatically.  The  two  primary  areas  of  Interest  have  been  the 
flow  downstream  of  a  shock  wave  (l.e.  within  a  shock  layer)  and  Internal 
eiq>andlng  flows.  The  former  Is  of  practical  Interest  when  designing  high 
speed  aircraft  or  reentry  vehicles  and  the  latter  is  of  Interest  mainly  In 
propulsive  devices.  Nonequilibrium,  expanding,  quasl*one  dimensional 
streamtube  flow  is  the  focus  of  this  investigation. 

A  basic  phenomenon  associated  with  expanding  compressible  flows  at 
high  temperature  and  low  pressures  is  the  freezing  or  trapping  of  energy  In 
chemical  or  Internal  energy  modes.  Hence,  energy  is  withheld  from  the 
system.  Departure  from  equilibrium  results  with  the  frozen  state  extreme 
clearly  containing  an  excess  of  energy  associated  with  Individual  particles 
over  that  which  It  would  have  in  an  equilibrium  state  at  the  local 
temperature  and  pressure.  This  Implies  a  further  reduction  In  temperature 
which  in  turn  effects  the  nonequilibrium  rate.  Since  convergent-divergent 
nozzles  or  streamtubes  effectively  convert  internal  energy  into  kinetic 
energy,  which  in  propulsion  applications  directly  translates  into  thrust, 
nonequilibrium  effects  may  degrade  nozzle  performance  appreciably. 

Some  Insight  into  the  trends  for  convergent-divergent  nozzle  behavior 
follows  on  examining  , the  classical  quasl-one  dimensional,  Inviscld, 
compressible  equations  such  as  derived  In  reference  [!}.  For  supercritical 
(i.e.  a  transition  from  subsonic  to  supersonic)  flow,  the  critical  section 


is  at  the  geometric  minimum  of  the  one  dimensional  nozzle 


Downstream  of 
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the  critical  point  both  temperature  and  density  decrease  rapidly  with 
increase  in  gas  velocity.  Typical  cooling  rates  are  of  the  order  10^  to 
10^  R/sec  [2] •  Since  changes  in  temperature  affect  the  equilibrium  energy 
distribution  curve  (i.e.  Maxwell-Boltzmann  distribution)  and  the  physical 
mechanisms  to  redistribute  the  energy  occur  at  a  finite  rate,  noticable 
departures  from  physical  equilibrium  are  possible. 

However,  nonequilibrium  processes  are  nonlsentropic,  hence  the  perfect 
gas  descriptions  for  inviscid  compressible  flow  (e.g.,  Liepmann  and  Roshko 
[1] }  require  modification.  For  nonlsentropic  processes  the  local 
thermodynamic  state  is  not  only  a  function  of  streamtube  area,  but  also  the 
upstream  area  distribution,  i.e.  the  history  of  the  fluid  flow  upstream  of 
any  point  in  the  streamtube.  The  specific  nonequilibrium  effects  present 
in  a  particular  flow  field  depend  on  the  gas  cotiq>onents  and  the  temperature 
and  density  of  the  flow.  This  follows  since  all  such  physical  phenomena 
need  particle  interactions  to  maintain  equilibrium,  and  temperature  and 
density  levels  indicate  the  energy  and  number  of  such  interactions. 
Representative  types  of  possible  nonequilibrium  physical  phenomena  to  be 
found  in  an  expanding  flow  are:  translational,  rotational,  vibrational, 
chemical,  and  radiation.  For  any  energy  mode  to  maintain  local 
equilibrium,  the  characteristic  time  for  readjustment  by  collisions  must  be 
negligible  compared  with  the  time  required  for  the  fluid  to  experience  a 
significant  change  in  local  conditions  (characteristic  flow  time).  In 
general,  the  adjustment  of  translational  and  rotational  modes  requires 
relatively  few  collisions  to  maintain  equilibrium  as  compared  to  vibration 
or  a  chemical  reaction.  Hence,  translation  and  rotation  have  short 
characteristic  times  for  readjustment  by  collisions  and  thus  are  generally 
assumed  to  be  in  local  equilibrium.  It  is  only  in  the  presence  of  large 
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flow  field  gradients  that  the  characteristic  readjustment  time  becomes  of 
the  same  order  as  the  characteristic  flow  time.  This  paper  emphasizes 
nonequilibrium  due  to  chemical  reactions,  assuming  all  other  modes  to  be  in 
local  equilibrium.  However,  the  methodology  is  applicable  to  a  general 
nonequilbrium  process  which  may  be  modeled  and  evaluated  using  the  same 
numerical  procedures. 

Since  chemical  reactions  occur  as  a  direct  result  of  molecular 
collisions,  it  is  apparent  that  as  the  temperature  and  density  decrease  in 
an  expanding  flow  the  number  of  interactions  decrease  and  energy  may 
effectively  become  frozen  due  to  ’’incomplete"  chemical  reaction.  If  we 
consider  the  simple  case  of  an  elementary  reaction  involving  diatomic 
dissociation,  the  atomic  species  concentration  (oc)  is  a  measure  of  the 
energy  required  to  break  the  diatomic  bonds.  Such  "tied  up"  energy  is  then 
unavailable  for  conversion  into  kinetic  energy.  In  a  flow  process  this 
manifests  Itself  by  several  distinct  regions  which  are  associated  with  the 
freezing  phenomenon.  In  subsonic  regions,  the  characteristic  flow  time  is 
large  since  the  temperature  and  density  are  high,  and  the  atomic 
concentration  closely  approximates  equilibrium  distribution.  As  a  flow 
accelerates  in  a  subsonic  section,  the  concentration  departs  from 
equilibrium  but  still  maintains  near  equilibrium  conditions.  The  flow  may 
then  continue  to  accelerate  with  corresponding  decreases  in  temperature  and 
pressure.  At  some  point  in  the  streamtube  such  decreases  imply  that  the 
characteristic  reaction  time  becomes  much  larger  than  the  characteristic 
flow  time,  since  sufficient  particle  collisions  to  maintain  equilibrium  are 
no  longer  present.  When  this  occurs,  the  concentration  quite  rapidly 
departs  from  equilibrium  and  effectively  "freezes  out",  essentially 
withholding  energy  from  the  flow. 


From  a  practical  standpoint^  a  nozzles,  for  example,  can  be  designed 
80  as  to  minimize  aonequlllbrium  effects*  There  are  basically  four  major 
parameters  of  Interest  here  that  control  nozzle  performance  and  they  will 
be  discussed  In  some  detail  later.  Briefly,  they  are  nozzle  scale, 
reaction  rate  constant,  stagnation  density,  and  stagnation  temperature. 
Previous  studies  Indicate  [3]  that  order  of  magnitude  changes  In  a 
Qondlmenslonal  rate  parameter  ($)  are  needed  to  significantly  alter  the 
nonequillbrluffl  concentration.  Increasing  the  stagnation  pressure  (po)  does 
have  an  effect  on  both  the  Initial  dissociation  and  the  recomalnatlon  rate. 
However,  for  propulsion  applications  there  will  be  a  weight  penalty  for 
Increases  In  the  stagnation  condition.  For  this  reason  predicting  thrust 
loss  due  to  nonequilibrium  effects  Is  important  to  the  performance  of 
propulsive  devices. 

Since  the  nonequilbirum  system  of  equations  is  nonisentropic, 
virtually  all  problems  of  practical  Interest  are  too  complicated  to  solve 
analytically.  Therefore,  various  schemes  have  been  proposed  to  numerically 
integrate  the  system  of  partial  differential  equations.  The  numerical 
techniques  have  primarily  been  related  to  the  difficulties  with:  1) 
boundary  conditions  and,  2)  near  equilibrium  situations.  One  of  the 
earliest  methods  [3]  employed  a  standard  numerical  integration  technique 
(Runge-Kutta)  applied  to  a  system  of  steady  Euler  equations.  For  an 
initial  valve  type  problems  (e.g.,  a  supersonic  inlet)  the  Inlet  conditions 
were  specified  and  a  solution  computed  by  a  space  marching  method. 
However,  if  one  boundary  is  upstream  of  the  critical  station,  use  of  a 
steady  Euler  equation  description  implies  a  two  point  boundary  valve 
problem  involving  the  critical  mass  flow  (m^,)  and  the  inlet  boundary 
conditions.  The  critical  mass  flow  is  determined  at  the  station  where  the 
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local  flow  velocity  is  equal  to  the  local  frozen  speed  of  sound.  Unlike 
the  isentropic  case,  (111^)  cannot  be  determined  a  priori  because  the  flow  at 
each  point  depends  on  the  entire  upstream  history.  In  addition,  for 
nonequilibrium  flow  the  critical  point  is  downstream  of  the  geometric 
minimum.  Therefore,  the  critical  mass  flow  rate  and  location  must  be 
determined  during  the  integration  procedure.  This  problem  is  illustrated 
by  the  following  relation  [2]: 


( 1)  *  du  =  dA  “ 

Ca.  vA  "STx  P 

'  V  ' 

which  is  the  nonequilibrium  extension  of  the  classical  expression  [I]: 


dA 

for  isentropic  flow.  For  —  *0  (geometric  minimum)  the  local  velocity 

dx 

equals  the  local  sound  speed  (ap  or  a^)  for  isentropic  supercritical  flow. 

However,  equation  (1)  Indicates  the  critical  point  can  occur  downstream  of 
ha  da 

the  throat  since  ( - )( — )  is  generally  positive.  This  is  true  because 

php  dx 

hp  is  negative,  h^  is  positive  (except  at  very  high  temperatures),  and  in 

da 

an  accelerating  flow  —  is  negative  due  to  decreasing  temperature*  At 

d  X 

such  a  downstream  critical  station  the  velocity  is  equal  to  the  local 
frozen  sound  speed  which  is  unknown  Initially* 

Several  methods  have  been  devised  for  the  two  point  boundary  valve 
problem  when  using  steady  state  equations.  A  simple  procedure  consists  of 
essentially  several  guesses  for  the  mass  flow  rate,  using  the  limiting 
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frozen  and  equilibrium  mass  flow  rates  as  guides  for  upper  and  lower 
bounds*  Successive  Iteration  were  computed  until  supercritical  flow 
downstream  of  the  throat  was  calculated  (Bray  [3],  Hall  and  Russo  [4])*  A 
second  technique  involved  Integration  to  just  upstream  of  the  critical 
singularity  with  a  guess  for  m^*  The  solution  was  matched  with  a  series 
expansion  about  the  singular  point  to  obtain  a  corrected  valve  for  m^,  and 
the  procedure  was  repeated  until  convergence  resulted  [3] •  An  Inverse 
method  proposed  by  Eschenroeder  et  al  [3]  eliminated  the  necessity  for 

Iteration  by  specifying  the  equilibrium  density  (p*)  Instead  of  the  area 

* 

distribution  In  the  transonic  nozzle  section*  Equilibrium  mass  flow  (m^) 
was  assumed  and  the  area  ratio  (A(x))  was  determined  as  part  of  the 
solution*  Integration  was  continued  downstream  of  the  critical  point,  then 
direct  integration  methods  were  used  for  the  remainder  of  the  computational 
domain*  Still  another  technique  (Bray  [3])  suggested  patching  together  an 
upstream  equilibrium  solution  to  a  nonequilibrium  flow  downstream  of  the 
critical  point* 

A  major  difficulty  In  the  Integration  procedure  Is  that  typically  the 
upstream  boundary  point  Includes  a  subsonic,  near  equilibrium  region*  The 
departure  from  close  to  equllbrlum  states  exhlbltes  a  singular  perturbation 
behavior*  This  results  from  the  fact  that  the  nonequilibrium  parameter 
(e*g*  concentration)  alters  Its  state  In  proportion  differences  from 
equilibrium  and  universally  proportional  to  a  characteristic  time;  In  the 
near  equilibrium  region  this  Is  essentially  0/0,  l*e*  closely 
indeterminate*  If  a  ordinary  explicit  space  marching  scheme  Is  used, 
extremely  small  Increments  must  be  taken  to  maintain  numerical  stability 
and  is  somewhat  Impractical*  In  reference  [5],  Lomax  and  Bailey  provide  an 
excellent  summary  of  numerical  integ. itlon  methods  applied  to  the  singular 
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perturbation  problem.  One  may  use  small  perturbation  theory  in  the  near 
equilibrium  region  to  increase  step  size  [6].  Anothe^^  technique  developed 
by  Treanor  [7]  uses  a  linear  approximation  to  the  rate  equation.  The 
linear  equation  is  integrated  explicitly  with  the  results  being  expressed 
In  the  form  of  the  Runge-*Kutta  finite  difference  equation  plus  a  correction 
term.  Step  size  can  also  be  increased  by  converting  the  differential 
equations  to  integral  equations.  Bray  [3]  avoids  the  problem  by  using  his 
patched  equilibrium  solution  technique.  Finally  implicit  numerical 
techniques  are  suggested  in  reference  [5],  and  a  review  of  implicit  methods 
as  applied  to  reacting  systems  is  given  in  reference  [8] .  Implicit  schemes 
do  remove  the  step  size  constraint  of  explicit  schemes  but  are  less 
accurate  and  require  more  computer  time  per  iteration. 

In  recent  years  consideration  has  been  give  to  circumventing  the  two 
point  boundairy  valve  problem  by  means  of  an  unsteady  description.  The 
system  of  equations  is  then  hyperbolic  in  time  throughout  the  entire  domain 
of  Integration  in  contrast  to  the  elliptic  subsonic  region  when  steady. 
The  technique  has  been  applied  to  ideal  compressible  [9]  and  nonequilibrium 
[10]  flows  using  explicit  integration  schemes  (e.g.,  Lax-Wendroff  [11]  or 
MacCormack  [12]).  Unsteady  methods  allow  integration  in  time  to  a  steady 
state  coiiq>atlble  with  imposed  steady  state  boundary  conditions.  The 
critical  mass  flow  and  critical  section  location  are  then  found  in  the 
course  of  the  numerical  evaluation.  However,  the  unsteady  scheme  still 
must  converge  to  a  stea4y  state  solution  in  all  near  equilibritim  regions, 
and  this  limits  step  size  in  much  the  same  way  as  when  using  steady  Euler 
equations. 

It  is  the  present  intent  to  take  advantage  of  both  time  dependent 
Euler  equations,  which  eliminate  the  two  point  boundary  valve  problem,  and 


Implicit  schemes  which  ellmlaate  the  stability  time  step  constraint*  In 
Chapter  2,  the  physical  model  for  nonequilibrium,  compressible  flow  In  a 
streamtube  will  be  described,  including  assumptions  and  model  equations* 
Boundary  conditions,  thermodynamics,  and  rate  chemistry  are  also  discussed* 
In  Chapter  3,  the  suggested  implicit  numerical  scheme  is  described  and 
numerical  treatment  of  the  boundaries  is  discussed*  Chapter  4  covers 
numerical  solutions  that  illustrate  the  capability  of  the  implicit  scheme 
when  applied  to  nonequilibrium  flow*  Finally,  Chapter  5  discusses  possible 
iiq>rovements,  extensions  to  higher  dimensions  and  larger  systems  of 


reactions 


Chapter  2 

ANALYTICAL  AND  PHYSICAL  MODEL  AND  GOVERNING  EQUATIONS 


The  modeling  is  that  of  a  nonequilibrium  flow  In  a  quasl-one 
dimensional  streamtube.  If  all  processes  were  Isentroplc,  algebraic 
equations  would  describe  the  thermodynamic  changes  of  state  [13] . 
Thermodynamics  then  constrains  the  inter-related  state  variables  but 
Implies  a  one  to  one  correspondence  between  physical  and  dynamic  properties 
calculated  from  the  equations  of  motion.  For  nonequilibrium  flow  the 
thermodynamic  behavior  depends  on  the  detailed  flow  history.  As  mentioned, 
the  unsteady  Euler  equations  will  be  used  to  model  the  nonequilibrium  flow. 
They  comprise  a  hyperbolic  system  of  partial  differential  equations  which 
describe  an  inviscid,  adiabatic,  compressible  flow*  The  quasl-one 
dimensional  system  is  then: 


mass; 
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In  addition  to  the  standard  (homogeneous  medium)  conservation 
equations,  there  is  a  single  conservation  of  species  constraint.  The 
physical  modeling  required  by  this  phenomenon  Is  presented  in  the  following 


section 
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2.1  Chemistry /Thermodynamics 

The  species  equation^  In  fact.  Is  representative  of  a  number  of 
possible  nonequilibrium  effects:  for  example,  Internal  modes  such  as 
vibrational,  rotational,  or  electronic,  and  chemical  modes  such  as 
dissociation,  ignition.  Here  we  will  assume  that  chemical  dissociation  is 
representative  and  the  sole  nonequilibrium  flow  phenomenon.  This  implies 
that  all  other  energy  modes  are  either  in  equilibrium  at  the  local  flow 
temperature  (T)  or  maintain  a  fixed  sub*-state,  i.e.  frozen.  The  species 
concentration  effects  the  system  both  implicitly  and  explicitly.  That  is, 
pressure  and  internal  energy  are  now  functions  of  concentration  as  veil  as 
two  other  thermodynamic  properties,  such  as  density  and  temperature. 

A  relatively  simple  binary  gas  suffices  to  illustrate  the 
nonequilibrium  effects.  An  approximate  but  realistic  gas  model  described 
by  Lighthill/Freeman  [13],  [14]  as  an  '‘ideal  dissociating  gas,"  and 
contains  the  major  feature  of  a  binary  reacting  system  while  simultaneously 
eliminating  complex  but  non-enlightening  terms .  Ligh thill  simplified 
several  contributions  to  the  partition  functions  for  a  diatomic  gas.  A 
brief  discussion  of  the  basis  is  described  below.  Consider  a  diatomic 
dissociating  gas  to  be  composed  of  particles  (A)  which  include 
translational  and  electronic  internal  energy  modes,  and  a  combined  particle 
(A2}  which,  in  addition  to  the  above,  has  vibrational  and  rotational 
degrees  of  freedom.  The  factorized  partition  functions  of  each  are  of  the 
form 


The  translational  partition  function  for  (A)  atoms  or  (A2)  molecules  is 
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where  m  is  the 
function  Is 

(0 


\Va 

=  V  A  Tj 

mass  of  the  specific  particle* 


The  rotational  partition 


and  the  vibrational  partition  function  is  given  by 


I 


j _ _ 


The  electronic  partition  function  is  described  by  the  general  equation 


AT 


where  the  are  degeneracy  factors  for  the  energy  levels  Specific 
values  for  the  constants  for  a  particular  species  are  obtained  by 
spectroscopic  analysis  of  transition  between  energy  levels*  Use  of 
equations  (5)  -  (8)  in  the  mass  action  equation  for  a  symmetric  diatomic 
gas  leads  to: 


where  (a*)  indicates  equilibrium  concentration*  Lighthill  noted  that  for 
1000  <  T  <  7000^K  the  contents  of  the  brackets  in  equation  (9)  is 
essentially  constant^  and  he  suggested  that  it  be  taken  to  be  a  constant, 
characteristic,  dissociation  density  The  simplified  expression  is 


then 


The  thermal  equation  of  state  is  unaffected  by  the  constant  p^j 
assumption  since  pressure  depends  only  on  derivatives  of  logarithms  of  the 
partition  function  with  respect  to  volume.  Assuming  both  A  and  A2  to  be 
perfect  gases,  the  thermal  equation  of  state  in  equilibrium  is 

(„•)  ^  -  c  1  4  o<*) 

The  expression  (1  +  a*)  is  essentially  the  equivalent  of  a  mixture  gas 
constant* 

Since  the  caloric  equation  is  a  function  of  temperature  it  is  affected 
by  the  Ligh thill  approximation.  The  general  expression  for  equilibrium  is 
[13]: 


0^)  C.  =  f  JnQ  +  ] 

,  nL  J 

Using  the  simplified  mass  action  relation  (10),  this  reduces  to  simply: 


In  which  contributions  due  to  electronic  energy  have  been  assumed  to  be 
negligible  for  the  temperature  range  considered.  The  internal  energies  of 
individual  constituents  of  the  mixture  are: 


Neglecting  electronic  energy  and  on  assuming  that  the  molecular  vibration 


mode  la  half  excited: 
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Since  e  ■  +  (l*a)e^  Eq.  (13)  then  follows.  Note  that  in  this  model 

^  4 

the  ratio  of  specific  heats  for  a  non-dissociated  gas  (a  «  0)  is  -  ,  which 

7  ^ 

is  less  than  the  true  -  for  diatomic  molecules.  From  Eqs.  (11)  and  (13) 

the  equilibrium  enthalpy  is: 


The  preceedlng  discussion  of  equilibrium  chemistry  provides  a  basis 
for  the  nonequilibrium  model.  The  thermal  (11)  and  caloric  (13)  equations 
of  state  are  of  the  same  form  for  nonequilibrium  and  require  only 
replacement  of  a*  by  a.  A  finite  rate  associated  with  concentration 
changes  requires  replacement  of  the  algebraic  mass  action  equation  by  a 
differential  description  equation  (3d).  The  right  hand  side  of  this 
equation  represents  the  production  or  recombination  of  atoms.  For  a  binary 
diatomic  dissociating  gas 


(it)  wCp,T,«) 


Ci-iV34 


Freeman  [13]  noted  that  for  equilibrium 


which  is  identical  to  the  Ideal  model  equation  (10)  If  Is  replaced  by 


He  also  replaced  the  first  parenthetical  factor  on  the  right 


hand  side  of  equation  (16)  by  a  single  Arrhenius  rate  constant. 


\  jVt 


Therefore  the  source  term  becomes 


(It) 


0 -  -  J2-  ^ 
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2 m3  Equilibrium  Flow 

A  limiting  case  to  be  considered  Is  that  of  Idealized  "Identical** 
equilibrium.  In  this  limit  the  differential  species  conservation  equation 
may  be  replaced  by  the  equivalent  law  of  mass  action  (10).  The  order  of 
the  equation  system  decreases  by  one  and  local  equilibrium  Is  enforced 
Instantaneously  and  continuously.  Thus,  the  state  vector  and  coefficient 
matrices  of  Eq.  (19)  are  now 


V  = 


4  ' 


where: 
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Actually,  for  the  Isentroplc  frozen  or  equilibrium  limits,  the  entire 
system  can  be  reduced  to  algebraic  forms.  However,  for  couslstency  such 
limiting  cases  have  been  Integrated  here  using  the  same  unsteady  Euler 
equation  description. 
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2«^  Time  Scales 


When  the  flow  is  neither  frozen  (a  *  constant)  nor  in  equilibrium, 
there  are  two  time  scales  present  in  the  problem.  Multiple  time  scales  can 
lead  to  integration  difficulties  (see  Chapter  3).  Such  problems  are 
strictly  numeric  and  result  from  the  presence  of  widely  separated 
eigenvalues  in  the  discrete  equations.  The  two  scales  are  a  flow  time 
(TfiQ^)  and  a  chemical  relaxation  time  (Tci^em)*  flow  time 


T. 


L/, 


♦  Us,  ^ 

is  a  measure  of  the  physical  scale  of  the  nozzle  or  residence  time  of  the 
fluid  as  was  suggested  in  Eq.  (20).  The  chemical  relaxation  time 


Us) 


is  a  measure  of  how  rapidly  the  fluid  will  relax  to  an  equilibrium  state. 
The  rate  parameter  4  in  Eq.  (21)  is  simply  the  ratio  of  the  time  scales, 
'^flow/’^chem* 

This  process  can  be  understood  in  terms  of  a  linearized  model  for  a 
rate  equation;  say 

in  which  T  is  a  modified  tQ^em* 

In  general  a*  and  t  are  both  time  dependent,  but  for  constant  values 
the  decay  towards  equilibrium  exhibits  a  simple  exponential  behavior  (Fig. 
1): 

(j-i) 


-  c. 


Clearly  t  is  the  time  required  to  decrease  the  relative  difference  by  the 
factor  e""^,  and  smaller  t  Implies  a  faster  relaxation  process.  Of  course 
Che  assumed  constant  heat  bath  in  general,  is  not  true  along  a  streamtube. 
Thus  variable  temperature  and  density  imply  varying  and  T^hem  according 
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to  local  conditions,  and  hence  there  is  a  coupling  to  the  fluid  motion, 
^flow 

The  limit  -  «  1  is  the  frozen  case  in  which  there  is  no  time  to  react 

"^chem 

to  local  flow  changes;  alternatively,  this  is  equivalent  to  a  short  path 

“^f  low 

length.  The  opposite  extreme,  -  »  1,  approximates  local  equilibrium 

‘^chem 

and  corresponds  to  a  relatively  long  local  residence  time  for  a  particle. 
It  is  the  near  equilibrium  situation  for  which  neither  time  scale  can  be 
neglected,  but  the  chemical  relaxation  time  is  relatively  small,  that  lead 
to  numerical  difficulties. 


2.5  Boundary  and  Initial  Conditions 

The  boundary  conditions  are  extremely  important  and  require  special 
care  in  the  case  of  numerical  evaluations.  For  a  subsonic  inlet  two 
boundary  conditions  are  required,  and  are  taken  here  to  be  stagnation 
enthalpy  and  entropy.  Since  interest  is  in  the  steady  state,  only  steady 
state  boundary  conditions  are  applied.  The  total  enthalpy  is  a  function  of 
temperature,  velocity  and  concentration: 


^  (  V  ♦  't')  T  +  ©4  j  ^ 

The  differential  entropy  is  [15] 

-  Cl  ^  -  ^4^  In  fC 

'  T  ^  C 

The  limiting  expressions  for  frozen  and  equilibrium  extremes  are  [13] 

5  =  3  In  X.  •<  (^l  -  (i^  (j- 4'^ 


T  =  Const 
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Chapter  3 

NUMERICAL  ANALYSIS 


3.1  Stiff aess/Stabiltty  Analysis 

It  has  been  pointed  out  that  the  nonequlllbrlum  system  of  equation 
(19)  la  numerically  stiff.  Stiffness  Is  not  the  result  of  an  analytical 
technique  but  can  originate  from  the  discretization  of  a  continuum  equation 
set  when  multiple  time  scales  (e.g.  Tf^ow  '^chem)  present.  The 

phenomenon  can  be  demonstrated  by  writing  the  species  continuity  equation 
in  finite  difference  form  using  explicit  upwind  differencing  (u>0): 

uAt 

where  c  - - (Courant  number) 

At 


At 

d  -  — 

‘^chem 

The  last  term  In  equation  (31)  represents  the  time  il-‘earlz«t'  'orm  of  the 
chemical  rate  equation  [13].  The  stability  limits  of  this  explicit  finite 
difference  alogorlthm  can  be  investigated  using  Von  Neumann  stability 
analysis  [16].  Each  Fourier  component  of  the  solution  is  written  as: 

n 


V  e 


which  on  substituting  into  (31)  yields 


V 


r>4 1 


I  -  c  (  I  '  )  -  c/ 


Since  the  modulus  of  the  amplification  factor  must  not  exceed  unity,  l.e. 


G.I  s 
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therefore 
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(,35)  I  ¥  CoS  +  cl(c(-^^ill 

If  satisfied  for  all  0  and  0<l-cos0<2 

(3C)  ‘yc(ctci-i')  iO 

or 

(37“)  e  ^  I  -  d/2^ 

The  classical  upwind  constraint  follows  for  the  frozen  limit  (d-^0),  namely 

(30  ^  ^  ‘ 

However,  notice  that  positive  c  requires 

(Ji)  cl  i 

or 

(yo)  At  i  - 

^  . 

Therefore,  as  (higher  reaction  rates)  the  step  size  (At)  for 

stability  is  determined  by  d  rather  than  c.  Hence,  for  situation  in  which 
step  size  is  dictated  by  chemistry,  the  number  of  iterations  for 
convergence  to  the  steady  state  may  be  impractical  for  explicit  techniques* 
An  implicit  method  is  of  some  interest  in  order  to  improve  the  step  size 
and  therefore,  iteration  count. 

3.2  Explicit/Implicit  Integration  Scheme 

The  specific  technique  applied  here  to  the  nonequilibrium  equation  set 
is  that  developed  by  MacCormack  [17]*  The  scheme  is  an  Implicit  analog  of 
MacCormack's  explicit  algorithm  [12]  introduced  in  1969*  In  fact,  the  new 
method  uses  the  earlier  explicit  technique  as  the  first  of  a  two  stage 
procedure  and  is  referred  to  as  an  explicit/implicit  algorithm.  The 
explicit  first  stage  uses  the  original  predictor  -  corrector  concept  to 
compute  local  changes  in  the  dependent  variables  of  the  governing 


equations •  Then,  the  implicit  second  stage  distributes  those  changes 


globally  to  all  points  in  the  computational  domain,  thus  removing  the 
stability  constraint  associated  with  the  explicit  algorithm.  The  procedure 
is  second  order  accurate  in  time  and  space,  unconditionally  stable  in  time, 
and  block  bidiagonal  in  form. 

3.2.1  Model  Equation 

The  basic  theory  and  implementation  of  the  explicit/implicit  algorithm 
can  be  demonstrated  using  a  model  convection  equation: 


^  +  e  ^  ■=  o 

The  explicit  part  integrates  equation  (Al)  using  a  predictor/corrector 


method: 
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U:  +  u; 
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In  the  predictor  step,  a  known  field  at  time  t  *  nAt  is  used  to  provide  an 
initial  estimate  for  the  new  field  at  time  level  t  «  (n+I)At  by  using  one 
sided  spatial  difference  to  compute  the  dependent  variable  change,  Au,  for 
the  hyperbolic  equation.  The  estimates  are  used  in  the  corrector 
step  with  opposite  one  sided  differencing  to  compute  a  second  estimate  of 
the  change  The  new  prediction  is  then  the  previous  value  plus  an 
average  of  the  two  estimates.  This  can  be  clearly  demonstrated  by 
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rewriting  the  finite  difference  equations  (42)  as  an  Euler  predictor 


followed  by 

a  modified 

Euler 

corrector  [18] : 

rv 

(»3) 

Ui  s 

+  At 

These  are 

equivalent 

to 

MacCotmack  when  the  time  derivatives  of  a 

hyperbolic  equation  (u*)  are  computed  using  forward  spatial  differences  in 
the  predictor  and  backward  spatial  differences  in  the  corrector.  Hence,  by 
equation  (43)  the  change  in  the  dependent  variable  (u)  is  the  average  of 
the  changes  at  the  predictor  and  corrector  levels.  The  method  is  second 
order  accurate  and  stable  if  the  time  step  satisfies  the  usual  CFL 
condition: 

Since  C  in  £q.  (41)  is  the  slope  of  the  characteristic  in  the  physical 
domain,  this  restriction  simply  insists  that  the  computational  domain  of 
dependence  must  include  the  physical  domain  of  dependence  (figure  2). 

The  stability  constraint  is  removed  by  Including  both  current  (t  » 
nAt)  and  new  (t  *  (n+l)At)  time  levels  in  the  finite  difference 


approximation  of  equation  (41).  Along  a  negative  slope  (c<0),  forward 
spatial  differencing  is  appropriate  since  the  information  is  propagating 
from  larger  x  at  earlier  time.  Therefore  an  implicit  form  for  equation 


(41)  becomes 

ur'= 


-  <  “c'") 


The  constant  (a)  essentially  determines  the  "degree  of  implicitness"  in  the 
equation.  If  (a  *  0),  the  original  explicit  equation  is  recovered;  if 
(a  *  1),  equation  (44)  is  fully  implicit.  On  collecting  terms.  Equation 
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(44)  can  be  written  In  the  following  form: 


where 
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For  stabllltry,  it  can  be  shown  that 

X  >,  t 


Note  that  the  first  term  on  the  right  hand  side  of  equation  (45)  is  the 
explicit  approximation  (42)  to  the  model  equation.  Hence,  the  implicit 
extension  of  the  original  MacCormack  explicit  scheme  Is 


c: 


A  t  X  = 
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+  4 


31 


This  system  is  unconditionally  stable  for  Integration  in  time  if  X 
satisfies  the  constraint  equation  (48).  In  addition,  both  predictor  and 
corrector  procedures  Involve  bidiagonal  matrix  equations  so  that  only  a 
single  sweep  through  the  computational  domain  is  necessary  for  each  level. 
The  bidiagonal  structure  can  be  seen  from: 
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3.2«2  Euler  System 

Using  the  basic  algorithm,  the  scheme  is  applicable  to  the 
nonequilibrium  Euler  System  (19)  as  represented  by: 


+  A 


H  - 


The  explicit/implicit  algorithm  for  matrix  equation  (49)  is 
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The  matrix  |a|  Is  related  to  the  eignevalues  of  the  system  by 


Cfi) 


where 
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The  matrix  of  eigenvectors  (S^)  satisfies 

■  A  =  s;'  s. 
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where  Is  the  dlagoaal  matrix  of  eigenvalues  of  (A)*  (Appendix  A 
provides  a  derivation  of  the  eigenvalues  and  eigenvectors  for  both  the 
equilibrium  and  nonequilibrium  equations) •  For  regions  which  satisfy  the 
usual  one  dimensional^  compressible,  explicit  stability  criteria 

(S-j)  4  t  6  4  iC 

I  U|  t  tx 

all  will  be  less  than  or  equal  to  zero*  In  such  cases  the  are  set  to 
zero  (see  equation  (52))  and  the  set  of  difference  equations  (50)  reduce  to 
an  explicit  form*  When  all  X^  are  positive  the  implicit  equation  can  be 
expanded  to  examine  the  matrix  block  bidiagonal  form*  For  example,  the 
implicit  predictor  equation  (50a)  can  be  written  as 

(rv)  (t  *  Ai  MirW  =4^7,%  At  Ml^,  sqr, 

\  y  Ax. 

and  solved  by  a  single  mesh  sweep  in  the  decreasing  (i)  direction*  The 
corresponding  corrector  sweep  is  in  the  opposite  direction*  Appendix  A 
provides  a  detailed  description  of  the  block-bidiagonal  solver* 

This  explicit/implicit  method  has  several  significant  advantages  over 
other  implicit  methods*  One  major  advantage  is  the  block-bidiagonal 
structure*  Numerous  other  implicit  schemes  (e*g*  Beam  and  Warming  [19]) 
eiiq>loy  a  block-trldlagonal  solution  algorithm  which  tends  to  increase 
coiiq)uter  time  per  iteration*  In  addition  the  bidiagional  method  is 
straightforward  to  program  and  does  not  require  matrix  inversion. 

An  area  of  difficulty  is  the  inclusion  of  a  source  term  (H)  in  the 
iiiq>llclt  stage  of  integration*  For  an  ideal  gas,  the  quasi-one  dimensional 
source  term  is  not  time  dependent  and  therefore,  does  not  effect  stability. 
However,  the  nonequilibrium  source  term  contains  the  time  dependent 
reaction  rate  (w)*  (See  equation  (19))*  As  will  be  shown  in  Chapter  4, 
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the  Implicit  stage  of  integration  will  provide  sufficient  stability  for 
moderate  reaction  rates.  But,  sufficiently  high  reaction  rates  effect 
stability  since  w  is  only  evaluated  during  the  explicit  stage  of 
integration. 

3.3  Boundary  Conditions 

The  solution  procedure  for  the  block  bidiagonal  matrix  equation  (48) 
assumes  that  all  dependent  variables  are  known  at  the  boundaries.  For 
certain  types  of  boundaries  (such  as  supersonic  inlets  or  outlets)  the 
boundary  conditions  are  straightforward.  However,  for  supercritical  nozzle 
flow  the  boundary  conditions  require  special  attention. 

Characteristic  theory  provides  a  basis  for  the  construction  of  stable 
Implicit  boundary  conditions.  As  noted  in  reference  [20],  a  correct 
formulation  of  boundary  conditions  are  of  extreme  importance  as  they  are 
the  governing  elements  of  the  computation.  Porous  boundaries  in  hyperbolic 
problems  are  represented  schematically  in  Figure  3  for  one  dimensional 
ideal  gas  flow.  (Reacting  flows  have  one  additional  characteristic;  see 
Appendix  A).  At  a  subsonic  inlet,  two  positive  eigenvalues  (Xx  *  u,  X2  * 
u-h))  define  paths  along  with  information  is  transmitted  to  the  boundary 
from  outside  the  integration  domain.  Consequently,  two  boundary  conditions 
are  needed  [21].  Along  the  remaining  negative  characteristic,  information 
is  transmited  from  inside  the  domain;  forward  spatial  differences  then 
correctly  account  for  the  Influence  of  such  characteristics  at  the 
boundary • 

To  apply  these  characteristic  constraints  to  the  implicit  algorithm, 
the  following  time  linear  approach  was  used  [22]  •  The  boundary  conditions 

(«)  Di  (  u)  =  o 

Can  be  written  in  time  linear  form; 
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Recall  (Chapter  3)  that  the  boundary  conditions  for  a  subsonic  Inlet  are 
that  total  enthalpy  and  entropy  are  fixed  at  the  Inlet.  Equation  (49)  In 
characteristic  form  Is 


is?') 

where 

At  the  Inlet, 
are  replaced 
becomes 


pU^+  APU^t  PH-O 
A  =  P'*  ^  P 

the  eigenvectors  corresponding  to  the  positive  eigenvalues 
by  the  boundary  conditions  (56).  Therefore  equation  (57) 
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Equation  (19)  Is  a  modified  partial  differential  equation  to  be  solved  at 
the  boundary  using  the  explicit/lmpliclt  method  (50).  Equation  (50)  need 
not  be  modified  at  the  supersonic  outlet  since  there  eigenvalues  are 
positive.  Implying  that  all  Information  at  this  boundary  Is  transmitted 
outward  from  the  Interior  of  the  domain.  Thus,  the  boundary  condition 
procedure  Involves  the  use  of  forward  difference  at  the  Inlet  and  backward 
differences  at  the  outlet.  The  Implementation  of  this  procedure  Is 
represented  schematically  by  Figure  4.  Notice  that  the  mesh  sweep  was 
applied  In  alternating  directions  to  eliminate  oscillations. 


Chapter  4 

NUMERICAL  SOLUTIONS 


4.1  General  Description 

A  computer  program  has  been  developed  from  application  of  the  physical 
modeling  and  numerical  techniques  of  Chapters  2  and  3  to  a  quasi-one 
dimensional,  nonequilibrium,  supercritical,  streamtube  flow*  Code 
validation  and  case  studies  (operating  on  a  PDF  11/23  computer)  are 
presented  here  to  indicate  the  utility  of  the  code  and  some  advantages 
relative  to  earlier  algorithms* 

Memory  limitations  (28  Kilobytes)  required  a  sequentation  of  the 
domain  in  order  to  achieve  larger  area  ratios*  An  upstream  segment  (see 
Figure  5)  contained  the  streamtube  minimum  section  and  required  the 
subsonic  inlet  boundary  condition  discussed  in  Chapter  3*  Once  supersonic 
flow  velocity  was  achieved,  a  fixed  supersonic  boundary  condition  was  used 
at  the  inlet  of  any  additional  downstream  streamtube  sections*  This 
procedure  allowed  for  expansion  to  an  arbitrarily  large  area  ratio  without 
either  exceeding  the  memory  limitations  of  the  system  or  sacrificing 
accuracy  as  would  be  required  by  an  otherwise  decrease  in  the  number  of 
mesh  points  for  a  given  area  change. 

Initial  conditions  (t*0),  boundary  conditions  (inlet,  outlet),  and 
coiiq>utational  methods  were  explored  in  the  selection  of  example  cases  for 
numerical  evaluation.  Table  1  indicates  the  individual  parameters  and 
specifications  for  the  study,  but  only  some  combinations  of  the  overall 
possibilities  were  actually  used  as  a  basis  for  the  solutions.  A  majority 
of  cases  were  considered  with  the  larger  stagnation  species  concentration 
(referred  to  as  the  "hot  gas"  case).  However,  a  selected  set  of  results 
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for  the  lesser  (''cold  gas”)  species  concentration  are  also  presented.  The 
area  distribution  along  the  streamtube  was  chosen  to  be  parabolic,  the 
larger  area  ratio  level  refers  to  the  exit  maximum  of  any  extended 
supersonic  section  (Figure  5).  Initial  conditions  are  the  assumed 
distributions  of  dependent  physical  property  variables.  Typically,  a 
decreasing  (downstream)  linear  variation  was  employed  in  order  to  be 
consistant  with  the  trends  of  expanding  flow. 

4*2  Numerical  Validation 

It  is  informative  to  compare  any  results  with  the  corresponding 
expected  trends  from  analytical  considerations  such  as  described  in 
references  [1],  [13]  or  [23].  Some  Important  checks  for  consistency  are: 
constant  mass  flow  rate  along  the  streamtube,  achievement  of  sonic  velocity 
at  the  minimum  section  in  the  limiting  frozen  or  equilibrium  calculations, 
and  appropriate  shifting  of  the  sonic  point  to  a  downstream  location  for 
nonequl librium  flow*  Of  course,  accuracy  of  a  numerical  result  depends 
upon  the  number  of  mesh  points  in  the  computational  domain,  integration 
step  size,  and  the  enforced  convergence  criteria.  The  latter  was  taken 
here  to  be 

M  ^  It?'”"’-  u  "I  <  ^  w*') 

as  the  criteria  for  all  computed  results. 

For  frozen  flow,  computed  results  were  compared  with  tabulated, 
analytical  results,  e.g.  [23].  Table  2  summarizes  the  numerical  solution 
of  an  ideal  gas  flow  (y  *  1*4)  and  compares  the  solution  with  the  exact 
analytical  values.  The  important  conclusions  are  first,  that  the  accuracy 
of  the  code  is  demonstrated  by  excellent  agreement  between  computed  and 
analytical  solutions*  That  is,  the  maximum  error  for  any  variable  at  any 
CFL  number  is  3*02  with  the  average  et.  r  generally  below  0*52* 
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Second,  the  characteristic  boundary  conditions  discussed  in  Chapter  3  do 
sat isf act orally  "find”  the  correct  set  of  Inlet  variables  such  that 
supercritical  flow  Is  achieved.  Third,  the  computed  Mach  number  at  the 
throat  deviates  from  unity  by  a  maximum  of  0.3%.  Finally,  the  number  of 
iterations  decrease  with  increasing  CFL  (i.e.  increasing  At),  but  with  an 
expected  decrease  in  accuracy. 

Numerical  solutions  for  equilibrium  calculations  were  validated  by  a 
coi]q>arison  between  computed  results  and;  the  law  of  mass  action  (10)  for 
consistency,  an  equilbriiom  sonic  speed  at  the  throat,  and  for  constant  mass 
flow  rate  along  the  streamtube.  In  all  cases  where  mass  action  was  checked 
the  agreement  between  the  numerically  computed  concentration  and  an 
analytically  computed  concentration  from  equation  (10)  was  accurate  to  six 
significant  figures.  The  equilibrium  Mach  number  departed  from  unity  by  at 
most  0.15%  (for  a  CFL  number  equal  to  10).  As  In  the  case  of  frozen  flow, 
the  accuracy  of  the  numerical  solution  improved  with  an  increase  in  mesh 
points  or  a  decrease  in  CFL  criteria.  For  all  cases  the  steady  state  mass 
flow  rate  varied  by  at  most  0.15%  and  for  nonequilibrium  flow  the  flux  had 
the  same  accuracy  as  for  isentropic  flow  at  a  given  CFL  number. 

In  addition  to  mass  flow  consistency,  nonequilibrium  flow  was 
validated  by  verifying  that  the  sonic  point  moved  downstream  of  the  minimum 
section  as  the  rate  coefficient  (<()  increased  (Equation  1).  Also 
qualitative  behavior  was  compared  with  the  results  of  Reference  [3]«  For 
all  cases,  the  equilibrium  and  frozen  limit  data  banded  that  of 
nonequillbrium  flow.  For  example,  an  increase  In  rate  coefficient  causes 
the  species  concentration  to  depart  from  near  equilibrium  *  farther 
downstream  and  the  flow  temperature  to  increase.  It  appears  from  the  cited 
general  results  that  the  computer  code  does  output  results  which  are 
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consistent  with  the  physics  of  the  problem. 

4.3  Integration  Difficulties 

Consider  problems  related  to  numerical  instability  as  specifically  due 
to  the  explicitly  modeled  source  term.  The  species  equation  source  term, 
w,  was  evaluated  at  the  current  time  level  as  outlined  in  Equation  (50). 
Such  a  basic  for  estimation  of  the  local  change  In  concentration  due  to 
chemical  reactions  can  be  extremely  unaccurate  when  the  rate  coefficient 
becomes  large.  Figure  6  illustrates  the  relative  scale  of  the  reaction 
rate  compared  to  temporal  step  size.  As  the  reaction  coefficient 
Increases,  some  flow  regions  tend  to  remain  in  near  equilibrium,  as  is  made 
clear  by  the  slope  of  the  fast  reaction  rate  in  Figure  6.  Essentially,  in 
such  "fast**  cases  any  small  deviation  from  equilibrium  produces  a  large 
Initial  rate  to  bring  the  concentration  back  to  an  equilibrium  level.  This 
initial  rate  quickly  decreases  as  the  chemical  relaxation  follows  an 
exponential  decay.  If  the  integration  step  is  a  sufficiently  small 
fraction  of  the  relaxation  time  (tchem)  such  as  used  in  Reference  [10],  the 
initial  rate  (or  slope)  will  decrease  at  the  next  iteration  level  due  to 
then  updated  values  for  dependent  variables.  However,  if  the  time  step 
exceeds  the  characteristic  relaxation  time,  the  concentration  may  (and  most 
frequently  does)  overshoot  the  equilibrium  level,  thus  leading  to  numerical 
instability. 

For  all  nonequilibrium  explicit  cases,  the  basic  predictor/corrector 
scheme  was  unstable  unless  the  time  step  was  taken  as  some  fraction  of  the 
smallest  relaxation  time  in  the  mesh.  This  constraint  is  apparent  from  the 
stability  restriction  of  Equation  (40). 

Vrhen  computing  using  the  explicit /implicit  algorithm,  the  local 
concentration  can  easily  (and  frequently  does)  overshoot  the  local  current 
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equlllbirum  level  without  achieving  numerical  instability  for  rate 
coefficient  magnitudes  less  than  $  *  10^.  A  possible  explanation  is  that 
the  implicit  stage  of  integration  modifies  the  explicitly  computed  local 
evaluation  by  distributing  these  changes  globally,  thus  effectively  damping 
out  oscillations.  Nevertheless,  as  the  rate  coefficient  increases  such 
numerical  damping  is  insufficient  to  maintain  stability. 

Four  distinct  methods  were  employed  in  attempting  to  improve  numerical 
stability.  Two  methods  used  a  cut  off  criteria  to  prevent  the 
concentration  from  overshooting  an  equilibrium  level.  A  third  method  used 
a  point  by  point  "type”  splitting  calculation  for  equilibrium  and 
nonequilibrium  regions.  That  is,  either  equilibrium  or  nonequilibrium 
equation  were  locally  solved  depending  upon  a  type  splitting  criteria.  A 
fourth  technique  consisted  of  a  modification  of  an  analytical/numerical 
approximation  suggested  by  Bray  [3] •  A  brief  description  of  these 
algorithms  will  be  presented  in  chronological  order. 

The  "type  splitting”  (Method  I)  procedure  employed  a  crossover 
criteria  to  determine  which  equation  set  (i.e.  equilibrium  or 
nonequilibrium)  was  to  be  solved  at  each  mesh  point  (Figure  5).  If  the 
local  reaction  rate  at  a  given  grid  point  was  sufficiently  large  such  that 
the  change  in  concentration  exceeded  the  change  to  the  current,  local 
equilibrium  level,  identical  equilibrium  was  assumed  for  that  mesh  point 
and  time  step.  Therefore,  equilibrium  Equations  (22)  were  applied  at  such 
points  and  nonequilibrium  equations  (19)  at  all  other  points  in  the  mesh. 
This  procedure  met  with  some  limited  success.  With  this  technique  the 
explicit  stability  limit  was,  in  fact,  substantially  increased.  By 
eliminating  the  explicit  relaxation  time  constraint,  the  temporal  step 
limitation  was  determined  by  the  CFL  criteria.  This  algorithm  proved  to  be 


stable  for  reaction  rates  up  to  ♦  *  10^*  To  add  an  implicit  stage  (and 
hopefully  further  increase  the  system  stability  limit)  required  an 
assumption  be  made  for  changes  in  concentration  (Aa)  at  "equilibrium** 
points.  For  identical  equilibrium^  a  is  completely  specified  by  the  lav  of 
mass  action*  However,  for  this  method  any  given  mesh  point  in  the 
computational  domain  could  be  in  either  equilibrium  or  nonequilibrium  at 
different  time  levels.  This  required  that  a  concentration  change  be 
specified  regardless  of  the  point  type.  Applying  an  estimated  change  did 
not  Improve  stability.  In  fact,  this  method  did  not  yield  any  stable 
Implicit  results. 

A  second  method  (II)  employed  an  analytical  approximation  by  Bray  [3] 
to  essentially  remove  the  singular  perturbation  region.  Bray  assumed  that 
in  regions  of  the  flow  field  for  which  the  dissociation  rate  term  was 
significantly  larger  than  the  species  convection  term,  the  flow  was  "In 
equilibrium. **  That  is,  the  suggested  criterion  is: 

(‘0  f  "V  ( I  - 

) 

Bray  considers  a  lower  limit  of  20  as  the  level  at  which  nonequilbrlum 
should  be  appropriate.  A  starting  point  was  chosen  based  on  this 
criterion.  Figure  7  indicates  that  the  procedure  is  not  unstable,  but  does 
result  in  an  oscillatory  solution  near  the  starting  point.  It  was  decided 
that  the  oscillations  were  related  to  the  near  equilibrium  nature  of  the 
starting  point,  which  lead  to  overshoot  for  large  time  steps. 

The  remaining  two  algorithms  (III, IV)  employed  explicit  damping  of  the 
species  equation  in  order  to  improve  and  control  stability.  The  first  of 
these  (Method  III)  used  a  cut  off  criterion  based  on  the  local,  time 
dependent,  equilibrium  concentration.  Any  concentration  change  that 
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exceeded  the  change  to  the  current  local  equilibrium  level  was  restricted 
(cut  off)  so  as  to  proceed  at  some  fraction  of  the  maximum  possible  change 
(a*  "  procedure  produced  numerically  stable  results  for  large 

time  steps  (CFL  *  10)  but  did  produce  some  ’’chatter”  that  could  not  be 

damped.  The  results  of  Figure  8  indicate  the  oscillations  in  a  typical 
solution. 

The  final  algorithm  (Method  IV)  eliminated  the  chatter  (of  III)  by 
basing  a  cut  off  criterion  on  the  local  steady  state  equilibrium 
concentration  in  contrast  to  the  local,  transient  equilibrium.  This 
procedure,  in  effect,  provided  a  lower  bound  on  the  concentration 
distribution,  and  added  the  necessary  stability  to  permit  consideration  of 
nonequilibrium  flows  with  large  time  steps  (Figure  9).  The  next  section 
provides  results  that  indicate  the  significant  improvement  that  becomes 
available  for  nonequilibrium  flow  evaluation. 

4.4  Results 

Method  IV  algorithm  implementation  was  as  follows.  First,  the 
computation  of  both  frozen  and  equilibrium  limiting  cases  established  an 
upper  and  lower  bound  to  the  nonequlllblrum  distribution  behavior.  The 
finite  rate  chemistry  cases  were  then  solved  using  an  initial  (t  0) 
concentration  distribution  selected  to  be  greater  than  the  equilibirum 
lower  bound  over  the  entire  domain.  During  the  explicit  portion  of  the 
algorithm,  if  an  implied  concentration  change  was  larger  in  magnitude  than 
that  required  to  achieve  local  steady  state  equilibrium,  the  allowed  change 
was  "cut  off”  at  some  fraction  of  the  difference.  That  is 
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The  entire  domain  was  scanned  at  each  step  to  locate  the  most  downstream 
point  for  which  a  cross  over  of  the  equilibrium  "stability”  boundary  would 
result  if  not  constrained  by  the  cut  off  condition  (Equation  (63))*  The 
cut  off  rule  was  then  applied  at  that  point  and  all  proceeding  upstream 
points*  By  this  means*  an  extremely  inaccurate  estimate  of  the 
concentration  change  was  inhibited  in  the  near  equilibriiun  region.  This 
effectively  restricted  the  consequences  of  the  local  reaction  rate  so  that 
the  transient  concentration  distribution  remained  above  that  for  the 
equilibrium  limit.  Essentially*  imposing  the  lower  bound  provided  the 

additional  stability  for  the  explicit/implicit  algorithm  in  the  singular 

. 

perterbation  region.  That  is*  since  the  chemical  reaction  rate  (w)*  is 
time  dependent  and  only  evaluated  at  the  current  time  level  the  inaccuracy 
of  the  explicit  estimate  of  the  rate  is  compensated  by  using  a  known  lower 
bound  to  inhibit  overshoot  and  oscillations.  The  individual  solutions 
provide  some  evidence  of  the  utility  of  Method  IV. 

The  purpose  of  the  completed  sample  cases  was  to  explore  the  ability 
of  the  algorithm  to  produce  results  consistent  with  the  physics*  and  to 
demonstrate  the  relative  gain  or  penalty  in  coizq>utation  time.  Five  cases 
(i.e.*  parameter  sets)  were  considered  in  order  to  indicate  a  range  of 
applicability  and  determine  possible  regions  of  difficulty.  The  five  cases 
were  as  follows: 


0,103,5x103,10^105,« 

0,103,5x103,10'^,105,« 

0,103,5x103,10^,» 


0.9  0,10^,« 
5.0  0,10^,- 


0.40 
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In  all  cases,  a  segmented  solution  procedure  was  employed  using  57  points 
between  computational  boundaries.  The  first  segment,  contained  the  minimum 
section  which  expanded  the  fluid  to  supersonic  velocity  «  2.0).  The 
second,  completing  supersonic  section  expanded  the  flow  further  over 
2.0  <  y  <  lO.O.  The  large  area  ratio  was  introduced  so  as  to  demonstrate 
chemical  freezing  even  more  clearly.  The  convergence  criterion  (Equation 
(60))  was  e  »  10*^.  The  initial  conditions  for  all  computed  solutions  were 
the  same.  That  is,  for  the  first  segment  *  2.0)  a  decreasing  linear 
distribution  for  density,  velocity,  and  temperature  and  a  constant 
concentration  distribution.  When  comparing  iterations  to  convergence  for  a 
given  stagnation  condition,  the  same  initial  profile  was  always  employed. 
Table  3-7  summarize  iterations  and  cut  off  criteria  for  each  case.  When 
comparing  relative  near  time  note  that  an  explicit/implicit  iteration  is 
2.25  times  an  explicit  iteration. 

4.4.1  Case  1  -  Hot  Gas  Explicit  Solution  (CFL  »  0.9) 

Results  for  explicit  integration  (CFL  *  .9)  are  Illustrated  in  Figure 
10-14  and  iterations  and  cut  off  conditions  are  summarized  in  Table  3.  The 
concentration  profiles  of  Figure  11  demonstrate  the  freezing  phenomenon  as 
supersonic  flow  continues  to  expand.  In  addition,  this  figure  shows  the 
tendency  of  the  concentration  to  maintain  a  near  equilibrium  distribution 
for  larger  area  ratios  for  increasing  The  initial  spatial  departure 
from  equilibrium  is  more  easily  seen  in  Figure  11  which  is  an  exploded  view 
of  the  throat  region.  ,  The  temperature,  velocity,  and  density  follow 
expected  trends  for  nonequilibrium  flow.  That  is,  temperature  and  velocity 
increase  and  density  decreases  with  increasing  An  interesting 
consequence  of  applying  the  cut  off  algorithm  is  that  the  explicit 
integrate  stage  (generally  restricted  by  chemistry  dominated  flows)  may 
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then  be  integrated  using  only  a  CFL  stability  criterion*  More  will  be  said 
regarding  the  quantitative  improvements  after  the  implicit  results  are 
discussed. 

4.4.2  Case  2  -  Hot  Gas  Explicit/Implicit  Solution  (CFL  *  5.0) 

Some  numerical  results  obtained  for  this  explicit/implicit  case  are 
sumxoarized  in  Table  4.  The  graphed  distributions  (Figures  15-19)  of 
dependent  variables  (p,u,T,a)  indicates  excellent  agreement  with  the 
explicit  (CFL  ■*  0.9)  solution.  For  example  the  exit  concentrations  for 
both  ^  *  10^>10^  differ  by  less  than  O.IZ  for  either  Also,  from  a 
comparison  of  Figures  II  and  16  it  can  be  seen  that  the  concentrations  at 
the  Intermediate  y  *  2*0  are  Indeed  virtually  identical.  However,  the 
implicit  (CFL  »  5.0)  algorithm  converged  in  approximately  0.56  the  time 
required  for  the  explicit  (CFL  »  0.9)  method. 

4.4.3  Case  3  -  Hot  Gas  Explicit/Implicit  Solution  (CFL  »  10.0) 

When  the  CFL  number  was  increased  to  10.0  a  significant  decrease  in 
accuracy  was  noted  for  higher  rate  constants.  Table  5  lists  the  iteration 
and  cut  off  coefficient  (K)  for  the  results  plotted  in  Figures  20-24*  Note 
that  results  for  ^  >  10^  are  not  presented  for  CFL  10.0.  In  order  to 
obtain  a  sufficiently  accurate  solution  for  $  »  10^  a  cut  off  criterion  of 
0.2  was  necessary*  This  increased  the  number  of  iterations  to  300,  three 
times  the  number  for  CFL  *  5.0  (compare  Tables  4  and  5).  Figure  25 
Illustrates  the  substantial  inaccuracy  that  corresponds  to  a  ^  «  10^ 
evaluation  with  K  »  1.0  used  as  a  cut  off  basis.  The  plot  of  those  results 
(Inverted  triangles  in  Figure  25)  indicates  a  virtual  equilibrium 
concentration  distribution  out  to  y  «  2.0.  However,  both  Cases  I  and  2 
indicated  that  departure  from  equilibrium  occurred  at  a  substantially 
smaller  area  ratio  for  the  same  rate  parameter  (see  Figures  10  and  IS).  A 
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decrease  to  K  -  0.2  (larger  values  were  tried)  led  to  solutions  for  Cases  2 
and  3  which  were  nearly  Identical,  as  Indicated  by  the  upper  concentration 
distribution  In  Figure  23,  but  with  a  run  time  Increased  by  a  factor  of 
three. 


A  measure  of  the  error  for  ^  •  10^  and  the  range  of  CFL  numbers  Is 
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plotted  In  Figure  26.  The  error  Is  defined  as: 

X  .31 


U3) 


e,rrcr  si 


%  (  o  u  ") /r 

I;  (o'"-  u  'V/i  ]' 


The  CFL  *  5  case  converges  much  more  rapidly  than  the  CFL  *  0.9  case  (100 
uses  410  Iterations)  and  produces  virtually  identical  results.  Similarly, 
for  CFL  ■  10  with  a  K  »  1.0  cut  off  condition,  convergence  Is  quicker  than 
for  either  0.9  or  3.0,  but  as  discussed  above  (Figure  25)  the  results  are 
quite  inaccurate.  The  K  *  0.2  iteration  history  plot  (Figure  26)  indicates 
the  substantial  increase  In  Iteration  count  (from  60  to  300)  and  suggests 
that  CFL  *  10.0  would  be  less  efficient  than  an  explicit  case  for  ^  «  10^. 
That  is,  when  the  acceptable  error  level  Is  below  that  corresponding  to  the 
explicit  stability  limit  then  an  explicit  method  becomes  preferable  to  an 
implicit  method  (Figure  27).  The  conclusion  is  that  for  a  sufficiently 
large  CFL  number  the  inaccuracy  Introduced  by  the  implicit  method  produces 
a  sever  constraint  on  the  cut  off  condition.  Therefore,  at  sufficiently 
high  a  lower  CFL  number  is  favorable  for  a  more  accurate  solution  in 
fewer  iterations. 

4.4.4  Comparison  of  Computation  Time 

The  improvement  in  number  of  iterations  and  computer  time  can  be  seen 
from  Figures  28  and  29.  Results  for  Cases  1  and  2  are  compared  with  an 
explicit  calculation  employing  t^i^em  stability  constraint  Instead  of  the 
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cut  off  algorithm*  As  discussed,  for  reaction  dominated  flows  the  explicit 
time  step  is  determined  by  the  characteristic  relaxation  time*  (e*g*  for 
an  upwind  algorithm  equation  (40)  determines  the  step  size)*  In  practice, 
numerical  investigation  has  indicated  [10]  that  for  stability  in  a  reaction 
dominated  flow 

CW)  At  = 

where 

^  <  ' 

That  is,  for  stability  the  integration  must  proceed  at  some  fraction  of  the 
smallest  relaxation  time  in  the  field*  The  results  shown  in  Figure  28  for 
the  explicit  Integration  without  the  cut  off  criterion  were  computed  with  6 
■  0.3*  This  value  was  determined  by  decreasing  g  until  a  stable  solution 
was  obtained*  The  results  plotted  in  Figure  28  indicate  that  the  "Tchem** 
case  deviates  from  the  others  downstream  of  the  throat*  As  would  be 
eiq>ected  the  cut  off  algorithm  does  Introduce  solution  error  but  the 
results  differ  by  only  0*1%  at  y  «  2.0*  The  other  dependent  variable 
(p,u,T)  indicate  similar  results  with  the  maximum  error  at  any  mesh  point 
being  less  than  3*0% •  The  significant  improvement  in  computer  time 
realized  for  this  small  decrease  in  accuracy  is  summarized  below: 


Method 

Iterations 

Average  At 

Relative  Run  Time 

Explicit  ~  CFL  “  .9 

e  • 

0.3 

3470 

0.04 

1.0 

Explicit  -  CFL  -  .9 

K  - 

1.0 

430 

0.2 

0.15 

Implicit  -  CFL  *5.0 

K  - 

0.5 

100 

1.0 

0.08 

This  chart  illustrates  a  drastic  Immprovement  in  run  time  for  a  small 
decrease  in  accuracy*  The  iteration  history  plot  (Figure  29)  shows 
graphically  the  difference  in  convergence  rate* 
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4.4.5  Cases  4  and  5  -  Cold  Gas 

For  completeness,  results  for  two  CFL  numbers  with  lower  stagnation 
conditions  (corresponding  to  »  0.40)  are  illustrated  in  Figures  30-33 
for  CFL  ■  0.9  and  Figures  34-37  for  CFL  «  5.0  and  iteration  and  cut  off 
conditions  are  summarized  in  Tables  6  and  7.  The  graphs  and  tables 
indicate  quite  similar  qualitative  behavior  and  conclusions  with  regard  to 
the  numerical  algorithm  consideration. 

4.4.6  Boundary  Conditions 

Characteristic  boundary  conditions  were  presented  in  Chapter  2  and  the 
numerical  procedure  was  discussed  in  Chapter  3.  The  characteristic 
boundary  condition  provided  stable  inlet  condition  without  overspecifying 
the  problem,  which  otherwise  may  have  prevented  the  automatic  development 
of  a  proper  supercritical  flow.  As  discussed  in  the  section  on  code 
validation  (Chapter  4),  a  choked  flow  condition  was  achieved  using  this 
characteristic  treatment  of  boundaries.  Figures  38  and  39  indicate  the 
solution  of  the  inlet  velocity  to  a  steady  state  value  consistent  with 
critical  mass  flow  and  unit  Mach  number  at  the  minimum  section  for  an 
equlibrium  ■  «)  calculation*  Figure  39  is  an  enlargement  of  Figure  38 
for  iteration  count  less  than  100.  For  all  CFL  numbers  the  inlet  velocity 
overshoots  the  steady  state  value  with  subsequent  oscillation  to  a  final 
value. 
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Chapter  5 

CONCLUSIONS  AND  RECOMMENDATIONS 


Considerable  savings  In  computing  time  has  been  demonstrated  using  an 
expllclt/lmpllclt  algorithm  In  conjunction  with  a  cut  off  condition  on  the 
explicitly  computed  reaction  rate  (w).  For  results  presented  In  this 
paper,  the  expllclt/lmpllclt  algorithm  was  two  to  three  times  faster  than 
the  explicit  stage  along,  when  the  cut  off  condition  was  employed  for  both 
the  explicit  and  expllclt/lmpllclt  methods.  In  addition,  when  comparing 
this  algorithm  with  other  explicit  Integration  methods  for  stiff  equations 
such  as  presented  In  Reference  [10]  the  explicit  stage  was  seven  times  and 
the  expllclt/lmpllclt  algorithm  twelve  times  faster  for  the  particular  case 
Investigated  ($  «  10^).  MacCormack*s  expllclt/lmpllclt  method  [15]  was 
proven  to  work  for  a  system  of  stiff  equation  and  with  the  addition  of  the 
cut  off  algorithm  was  extended  to  cover  a  wider  range  of  nonequlllbltiim 
flow  (0  <  *  <  105). 

Characteristic  boundary  conditions  have  also  been  devised  for 
quasl-one  dimensional  superclrtical  flow.  These  boundary  conditions 
permitted  the  solution  of  an  unsteady  system  of  equation  (to  a  steady  state 
final  value)  without  ^  priori  specification  of  inlet  conditions  consistent 
with  supercritical  mass  flow.  Results  presented  Indicate  that  with  these 
characteristic  boundary  conditions  the  Mach  number  at  the  minimum  sections 
for  the  limiting  equilibrium  and  frozen  cases  was  accurate  to  0.152. 

An  extension  of  the  approach  may  be  fruitful  In  three  areas:  1) 
increased  number  of  dimension;  2)  higher  order  reaction  systems;  3) 
different  Integration  algorithms.  The  extension  of  the  current  method  to 
higher  dimensions  will  require  careful  consideration  of  the  stability 
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boundary  concept  when  applying  a  cut  off  condition*  It  may  be  possible  to 
simply  apply  the  cut  off  condition  in  an  analogous  manner.  That  is,  solve 
the  two  dimension  equilbrium  case  and  use  this  as  the  local  lower  bound* 
For  higher  order  reaction  systems  possible  concentration  stability 
boundaries  might  be  defined  by  considering  the  equilibrium  distribution  of 
each  species  with  all  other  species  frozen*  Finally,  for  improved 
Integration  techniques,  a  more  extensive  parametric  study  of  cut  off 
condition  with  reaction  rate  coefficient  and  CFL  number  must  be  performed 
to  provide  more  data  regarding  the  optimum  cut  off  criterion  combination  in 
terms  of  a  given  reaction  rate  coefficient  and  step  size.  In  addition,  as 
discussed  in  Chapter  3,  inclusion  of  the  source  term  in  the  Implicit  stage 
of  integration  should  be  Investigated  to  determine  if  numerical 
instabilities  due  to  explicit  evaluation  of  the  rate  term  can  be  eliminated 
without  losing  the  advantage  of  the  bidiagonal  structure  of  the 
explicit/implicit  algorithm* 

In  all,  substantial  improvements  in  computation  time  for  a  system  of 
stiff  equations  have  been  demonstrated.  Application  of  characteristic 
boundary  conditions  in  conjunction  with  an  unsteady  system  of  Euler 
equations  has  eliminated  the  two  point  boundary  value  problem.  The  cut  off 
algorithm  has  extended  the  range  of  applicability  of  the  MacCormack 
explicit/implicit  scheme.  The  success  of  this  numerical  method  provides  a 
basis  for  extension  to  higher  dimensions  and  more  complex  reacting 
systems. 


51 


REFERENCES 


[1]  Llepmann,  H.W*,  Roshko,  A.,  Elements  of  Gasdynamics*  John  Wlley»  New 
York,  1957. 

[2]  Hall,  J.G.,  Treanor,  C.E%,  ’’Nonequlllbrlum  Effects  In 
Supersonic-Nozzle  Flows.”  AGARDograph  124,  December  1967. 

[3]  Bray,  K.N.C.,  "Chemical  and  Vibrational  Nonequilibrium  in  Nozzle 
Flows."  Nonequilibrium  Flows  Part  II,  ed.  Peter  P.  Wegner,  Marcel 
Oekker,  Inc. ,  New  York,  1970. 

[4]  Hall,  J.G. ,  Russo,  A.L. ,  "Studies  of  Chemical  Nonequilibrium  in 
Hypersonic  Nozzle  Flows."  Cornell  Aeronautical  Laboratory,  Report 
AD-1118-A-6,  AFOSR  TN  59-1090,  1959. 

[5]  Lomax,  H.,  Bailey,  H. ,  "A  Critical  Analysis  of  Various  Numerical 
Integration  Methods  for  Computing  the  Flow  of  a  Gas  in  Chemical 
Nonequillbrium. "  NASA  TND-4109,  August  1967. 

[6]  Eschenroeder,  A.Q.,  Boyer,  D.W.,  Hall,  J.G.,  "Nonequilibrium 
Expansions  of  Air  with  Coupled  Chemical  Reactions."  The  Physics  of 
Fluids,  Volume  5,  Number  5,  May  1962. 

(7J  Treanor,  C.E.,  "A  Method  for  the  Numerical  Integration  of  Coupled 
First-Order  Differential  Equations  with  Greatly  Different  Time 
Constants."  Mathematics  of  Computation,  Volume  20,  1966. 

[8]  Kee,  R.J.,  Dwyer,  H.A.,  "Review  of  Stiffness  and  Implicit 

Finite-Difference  Methods  in  Combustion  Modeling." 

Progress  in  Astronautics  and  Aeronautics  Volume  76,  eds.  J.  Ray  Bowen 
et  al,  1979. 

[9]  Serra,  R.A*,  "The  Determination  of  Internal  Gas  Flows  by  a  Transient 
Numerical  Technique."  AIAA  Paper  71-45,  1971. 

[10]  Anderson,  J.D.,  "A  Time-Dependent  Analysis  for  Vibrational  and 
Chemical  Nonequilibrium  Nozzle  Flows."  AIAA  Journal,  Volume  8,  Number 
3,  March  1970. 

[11]  Lax,  P.D.,  Wendroff,  B.,  "Difference  Schemes  for  Hyperbolic  Equations 
with  High  Order  of  Accuracy."  Communications  on  Pure  and  Applied 
Mathematics.  Volume  17,  1964. 

[12]  MacCormack,  R.W. ,  "The  Effect  of  Viscosity  in  Hypervelocity  Impact 
Cratering  ."  AIAA  Paper  69-354,  May  1969. 

[13]  Vincent!,  W.G.,  Kruger,  C.H.,  Introduction  to  Physical  Cas  Dynamics. 
John  Wiley,  New  York,  1965. 

[14]  Llghthlll,  M.J.,  "Dynamics  of  a  Dissociating  Gas  Part  I;  Equilibrium 
Plow."  Journal  of  Fluid  Mechanics,  Volume  2,  1957. 


Journal 


[15]  Clarke,  J.F.,  ”The  Linearized  Flow  of  a  Dissociating  Gas.** 
of  Fluid  Dynanics*  Volume  7,  1960. 

[16]  Roache,  P-J.,  Computational  Fluid  Dynamics.  Hermosa  Publishers, 
Albuquerque,  1976. 

[17]  MacCormack,  R.W.,  "A  Numerical  Method  for  Solving  the  Equations  of 
Compressible  Viscous  Flow.”  AIAA  Journal,  Volume  20,  Number  9, 
September  1982. 

[18]  Kutler,  P. ,  Lomax,  H.,  ”Shock-Capturing,  Finite-Difference  Approach 
to  Supersonic  Flows.”  AIAA  Paper  71-99,  January  1971. 

[19]  Beam,  R.M.,  Warming,  R.F.,  ”An  Implicit  Factored  Scheme  for  the 
Compressible  Navier-Stokes  Equations."  AIAA  Third  Computational 
Fluid  Dynamics  Conference,  Albuquerque,  New  Mexico,  1977. 

[20]  Moretti,  G. ,  "Improtance  of  Boundary  Conditions  in  the  Numerical 
Treatment  of  Hyperbolic  Equations.”  The  Physics  of  Fluids, 
supplement  11,  1969. 

[21]  Gustafsson,  G. ,  Ollger,  J.,  "Stable  Boundary  Approximations  for  a 
Class  of  Time  Discretizations  of  u^  =  ADqu.”  OPPSALA  University 
Department  of  Computer  Science,  Report  Number  87,  September  1980. 

[22]  Chakravarthy ,  S.R.,  "Euler  Equations  -Implicit  Schemes  and  Implicit 

Boundary  Conditions."  AIAA  Paper  82-0228,  January  1982. 

[23]  Ames  Research  Staff,  "Equations,  Tables,  and  Charts  for  Compressible 
Flow."  NACA  Report  1135,  1953. 

[24]  Warming,  R.F.,  Beam,  R.M.,  Hyett,  B.J.,  "Diagonalization  and 
Simultaneous  Symnetrizatlon  of  the  Gas-Dynamic  Matrices . " 
Mathematics  of  Computation,  Volume  29,  Number  132,  October  1975. 


53 


Table  1 

INVESTIGATED  NUMERICAL  MODEL  PARAMETERS 


Method 

Explicit 

Explicit/Implicit 

Initial  Conditions 

Dependent  Variables 
(p,u,T,a) 

Configuration 

Maximum  area  ratio 

for  parabolic  distribution 


CFL  Number 
0.9 

5.0,  10.0 


*  Constant  stagnation  level 
Linear  between  spatial 
boundaries 

-  Equilibrium  distribution 

2.0351,  10.31 


Stagnation  Conditions 


Concentration 


0.4,  0.67 
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Location 

Inlet 


Throat 


Outlet 


Table  2 

MACCORMACK  EXPLICIT/IMPLICIT 
ALGORITHM  WITH  CHARACTERISTIC  BOUNDARY 
CONDITIONS 


Input 

data:  y  • 
No. 

^max  “ 
mesh  points 

2.0351, 
-  57,  e 

-  10-6 

CFL 

P. 

u 

I 

Max  X  Error 

Iterations 

Exact 

0.9564 

0.2973 

0.9823 

- 

- 

0.9 

0.9563 

0.2976 

0.9823 

0.1 

600 

5.0 

0.9565 

0.2975 

0.9824 

0.07 

120 

10.0 

0.9574 

0.2977 

0.9827 

0.13 

80 

Exact 

0.6339 

0.9128 

0.8333 

- 

- 

0.9 

0.6341 

0.9129 

0.8333 

0.03 

600 

5.0 

0.6337 

0.9114 

0.8330 

0.16 

120 

10.0 

0.6338 

0.9091 

0.8327 

0.40 

80 

Exact 

0.1806 

1.5743 

0.5044 

- 

- 

0.9 

0.1809 

1.5741 

0.5051 

0.17 

600 

5.0 

0.1817 

1.5730 

0.5081 

0.73 

120 

10.0 

0.1827 

1.5721 

0.5119 

3.0 

80 
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Table  3 

CASE  1  SUMMARY 

Oq  -  0.67,  CFL  -  0.9,  e  -  10"^ 


i 

2 

Y 

iinax 

Iterations 

K 

Iterations 

K 

0 

400 

- 

180 

- 

103 

410 

1.0 

180 

- 

5x103 

400 

1.0 

190 

- 

10^ 

410 

1.0 

200 

- 

105 

430 

I.O 

220 

1.0 

OP 

460 

Table  4 

210 

Oo  “ 

CASE  2  SUMMARY 
0.67,  CFL  -  5.0, 

e  -  10-^ 

i 

2 

Imaz. 

12. 

Iterations 

K 

Iterations 

K 

0 

90 

- 

60 

- 

103 

100 

1.0 

60 

- 

5x103 

90 

1.0 

60 

- 

10^ 

100 

1.0 

60 

- 

105 

100 

0.5 

70 

I.O 

100 


60 
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Table  7 

CASE  5  SUMMARY 
Oq  *  0.4,  CFL  ■  5.0,  e 


±  Jmax 

2 

^  Iterations  K 

< 

0  90  - 

10^  100  1.0 

«  100 


4 

J 

i 

4 

! 

4 

i 

i 


.  10"^ 

10 

Iterations 

60 

60 

70 


Characteristic  Boundary  Conditions 


amd  Cut-off  Condition 


concentration  profiles 
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RATIO 


CONCENTRATION  PROFILES 
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TEMPERPlTURE  PROFILES 


fiREA  RATIO 


VELOCITY  PROFILES 


OREft  RftTIO 


DENSITY  PROFILES 
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PREA  RATIO 


CONCENTRpiTION  PROFILES 


fiREfi  RATIO 


CONCENTRATION  PROFILES 


AREA  RATIO 


TEflPERATURE  PROFILES 


PREP  RPTIO 
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OREft  ROT 10 


DENSITY  PROFILES 


75 


0-X  8*0  3*0  3*0  0*0 


LiISN30 


PREP  RPTIO 


CONCENTRr^TION  PROFILES 


OREO  ROTIO 


TEflPERf^TURE  PROFILES 


3aniyy3dU3i 


PREP  RPTIO 


VELOCITY  PROFILES 


OREO  ROTIO 


DENSITY  PROFILES 


80 


0*T  8*0  3’0  2*0  0*0 


JLilSNSa 


PREP  RPTIO 


RATE-1. 0E* 
>--CFL-5.0 
^--CFL-10. 
•"•--CFL-10. 


81 


N0IiyyiN30N00 


AREA  RATIO 


ERROR  HISTORY 


ITER/^TIOHS 


83 


Nuiriber  of 


Ej^licit  Implicit 


Fig.  27.  Trends  in  Algorithm  Accuracy  and  Iteration  Count 
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Appendix  A 

EIGENVALUES,  EIGENVECTORS  AND  BOUNDARY  CONDITIONS 


Dlagonalzatlon  of  Euler  Equation 

As  noted  in  Chapter  3  the  system  of  nonequilibrium  unsteady  Euler 
equation  can  be  solved  using  a  block  bidlagonal  implicit  scheme.  In  order 
to  Implement  this  algorithm,  eigenvalues  and  eigenvectors  of  the  system 
must  be  computed.  This  section  will  cover  their  deviation  following 
procedures  similar  to  Reference  [24] • 

Nonequilibrium  Case 

For  the  nonequilibrium  case  the  system  of  equations  is: 


U- 

o 

o 


o 

u. 

o 

o 

o 

9 

u 

0 

T 

J 

.  J 

O 

SI 


16.  Id  t 
A 


or 

"^4  -t  >4  u.  +  H  =  O 

To  diagonalize  this  set  of  equations  the  eigenvalues  and  eigenvector  of  A 
must  be  found  from  the  relations 


(«)  Xj  A  =  Aj  X, 

or 
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C/fi)  <tCA-Xil')  =  0 

where  Xf  and  Xj^  are  the  row  eigenvector  and  eigenvalues  of  matrix  A. 
The  eigenvalues  are  computed  by; 

Us)  det  (a  '  At  l)  -  o 

which  becomes 

\  a))  (  A  -  1.  *  o 

(R+Cv) 

Here  a^  - - RT  which  Is  the  ratio  of  the  local  frozen  speed  of 

Yo^o^v 

sound  to  the  stagnation  speed  of  sound*  Hence  the  eigenvalues  are 

X,  =  lA-  a 

X^  *  K 

^  ^  —  tA.  ^  A 

\h  -  ^ 

The  eigenvalues  are  calculated  by  substituting  the  eigenvectors  Into 
Equation  (A *3) 


These  eigenvectors  form  the  rows  of  a  matrix  Now  the  matrix  A  can 

^e  diagonalized  by  the  similarity  tranform 

La>)  a  =  s.  a  s;' 


where  A  is  the  diagonal  matrix  of  eigenvectors 
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Equilibrium  Case 


For  the  equilibrium  system  of  equations: 


(515) 


Vx(2.  ->Tj 

fil  is 

^  *  =0 

-"‘0:?  Vx 


=  o 


tuK«rc  *^1  * 


By  the  same  procedure,  the  eigenvalues  of  A  are  here 


Kx  * 


3 


X,  -  u.  -  a 


M 


where 


Z 

Ci 


j  I  Cl  - <?‘0^  )  »(g^3o( 

f  *■  5Ca."<XT/<b;f  3 


vhich  is  the  ratio  of  the  local  equilibrium  sound  speed  to  the  stagnation 
sound  speed. 

The  matrix  of  eigenvalues  Is  then  computed  using  (A.5) 
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o 


e  r  R  -ni^T 
T  I  Rte„f><^. 


-  e 

_e.r  R  *t.,J  >^t1 

1*  [  K  +  1 


Boundary  Conditions 

The  Implicit  boundary  condition  at  the  subsonic  inlet  is  applied  as  in 
Reference  [22]  •  For  the  case  of  equilibrium  flow  there  are  two  positive 
and  one  negative  eigenvalue.  Therefore  two  boundary  conditions  can  be 
applied: 

S  «  constant 
hQ  *  constant 

For  an  ideal  dissociating  gas 


(4/c)  S  ^  "SinX  +  *<(  h(|- •<') 

-  ^  /  +  •^\  In  ^  f 

C/fl<7)  K.  =  R^|(V<-)T  ^ 

Appling  these  boundary  conditions  in  time  linear  form 

Vi  iv  =  o  Hit  -  o 
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Therefore, 


(fn) 

(»ii) 

(4^3) 


>e 

f  T 

Is-  = 

-2.  + 

>T 

T  f*  "iT 

Us.  = 

If 

K.  =  R,  C‘<+-.')  <•  (Tt©4^  U 


By  eliminating  the  eigenvectors  of  the  positive  eigenvalues  and  replacing 
them  with  the  boundary  conditions 


P 


1  -  f.  «.  ? 

J2- 

ft  *T 

T 

-5 

Cv  ^+lTf  G^V* 

-L  +  ^  hL 

?  T 

T  T  "^T 

• 

1  -aV.R.^ 

ftT 

1 

O  o 

o 

O  o 

o 
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Therefore,  the  system  of  equations  at  the  boundary  becomes 

t'WO  P,  ^  Pa  A  U.  H  -  O 

or  . 

C-U7)  -  -P.  -^H) 

The  ezpllclt/lmpllclt  algorith  Is  used  at  this  point  with  forward 
differences  only.  At  the  supersonic  outlet  all  eigenvalues  are  positive 
thus  no  special  treatment  Is  necessary. 

Algorithm  Implementation 

With  the  eigenvalues,  eigenvectors  and  boundary  conditions  specific, 
the  Implementation  of  the  Implicit  scheme  Is  as  follows: 

Recall  that  the  bldlagonal  matrix  equation  Is 


Rather  than  Invert  the  matrix 

The  following  matrix  algebra  may  be  performed 

or 

Th«  diagonal  matrix 


103 


proves  to  be  such  more  easily  integrated. 

The  cosKplete  sequence  is  as  follows:  (Reference  [17]) 

2)  X  '  W 

3)  0^  o»/  cw  l*f  ej 

4)  y  =  ( I  4  ^  o^y  X 

»  -  sr  y 

6)  Z  =  y 


7)  I  1  S'C'i’'" 
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Appendix  5 

COMPUTER  PROGRAM  INPUT  VARIABLE  NAMES 

The  following  Is  the  list  of  user  supplied  input  variable  names  for 
the  computer  program.  The  definition  of  some  geometric  variables  are 
illustrated  in  Figure  El*  Dimensions  are  given  where  appropriate. 


CAl 

Species  concentration  at  i  «  1 

CF 

Chemical  rate  constant 
(Reference  [13])  page  231) 

CFL 

Courant  number 

DELMAX 

Convergence  criteria  (Equation  (60) ) 

ETA 

Temperature  exponent 
(Reference  [13]  page  231) 

IBUG 

Debug  printout  if  ■  1 

ICHEM 

Nonequilibrium  if  *  2 

Equilibrium  if  *  3 

ICON 

Detail  printout  for  specified  mesh 
point  if  ■  1 

ICTOFF 

Apply  cut  off  criterion  if  *  1 

I6E0 

*  0  no  minimum  section 

*  1  minimum  section 

INLET 

*  0  characteritic  inlet  boundary 

conditions 

*  1  fixed  inlet  boundary  condition 


IPT  Specify  mesh  point  for  detailed  printout 

IREAD  «  1,  read  initial  conditions  from  unit  98 

ITAPE  •  1,  write  dependent  variables  of  final 

iteration  to  unit  99 
■  2,  write  error  history  to  unit  77 

ITMAX  maximum  number  of  iterations 

IWRITE  interval  for  writing  data  to  output  file 

NMAX  number  of  mesh  points 
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PO 

Stagnation  pressure  [Atm] 

RD 

Characteristic  dissociation 
density  [g/cm^] 

RBOl 

Initial  inlet  density 

RB02 

Linear  decay  coefficient 

THD 

Characteristic  dissociation 
temperature 

TBLX 

Reaction  rate  coefficient  ($) 

TO 

Stagnation  temperature  [^K] 

T1 

Initial  Inlet  temperature 

T2 

Linear  decay  coefficient 

VO 

Initial  Inlet  velocity 

VI 

Linear  growth  coefficient 

XKCHEM 

Cut  off  coefficient  (K) 

XL 

Distance  from  mlmlnum  section  to 
point  where  YB  Is  specified  (Figure  Al) 

XLT 

Total  length  If  IGEO  -  1 
length  to  exit  from  origin  if 

IGEO  -  0 

XKA2 

^folecular  weight  of  diatomic  species 

TB 

Half  height  at  X  *  XL 

■  -iW 


All  data  entered  is  formated.  Real  variables  are  read  using  FIO.O 
format,  except  for  TRLX  and  CF  which  are  read  in  E10.3  format.  Integer 


variables  (those  beginning  with  I  or  N)  are  read  with  15  format.  The  data 
entry  sequence  is  as  follows  with  each  line  representing  a  data  card  or 
line  entry: 


NMAX 

ITMAX 

IWRITE 

IBUG 

DELMAX 

XL 

XLT 

YB 

CFL 

RHOl 

RH02 

T1 

T2 

VO 

ICHEM 

TO 

PO 

XMA2 

RD 

THD 

ETA 

TRLX 

CF 

ICON 

IPT 

ICTOFF 

IGEO 

INLET 

IREAD 

ITAPE 

CAl 

XKCHEH 
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XLT 

(IGEO  =  0) 


Fig,  B  1.  Geometry  Input  Variable 
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Appendix  C 

PROGRAM  SUBROUTIt^E  DESCRIPTION 

The  following  is  a  list  of  subroutines  used  to  solve  the 

nonequilibrium  streamtube  flow  probvlem*  A  description  of  each  subroutine 

is  provided  below. 

MAIN 

-  executive  routine  to  manage  input,  output,  initialize  data, 
and  check  convergence. 

INTEG 

*  integration  subroutine.  Computes  explicit  and  Implicit  changes 
in  dependent  variable  Ul  and  applies  cut  off  criterion  if  ICTOFF 
»  1.  Called  from  PROGRAM  MAIN. 

EVBCT 

-  computes  matrix  of  eigenvectors  (SX)  and  its  inverse  (SXINV). 

Called  from  SUBROUTINE  INTEG. 

BNDRY 

*  calculates  the  characteristic  boundary  condition.  Called  from 
SUBROUTINE  INTEG. 

CHEM 

-  calculates  the  species  continuity  source  term  OMEGA.  Called 


from  PROGRAM  MAIN 
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Appendix  D 

COMPUTER  PROGRAM  LISTING 

The  following  pages  contain  the  complete  listing  of  the  nonequilibrium 
stteamtube  flow  program  described  in  this  paper. 
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VN  IV  V02.'‘'  Mon  24-Jc-.n-g3  PAGE  001 

PPCGRAh*  main 

COMMON  Ul(61»4)»iX{4.4),:xlNV(4i4)^DU(61^4)»U(61»4)^ 

1  A(61)iCiCrLirT*r'''*GAM0»II»NMA,x»au-»i;UGirLUXS(4)»ICALL 

2  » I  r  N  D 

3  »0M"GA(61)»TC(-"'‘(61)»Q(51)»ACH<=-M,cvrirTA»rCHGM»'’A2» 

4  PP»PHOO»rO-»TJ‘'',’rpLX»TO»CE(bl)OCrRl»DCDT10NP»IEORAT 

5  »ECnN(61)  »  ICTn"r->XKrHE’-'»lNLFT,N  ilN 

DIMENSKN  X  (6  1  )  ^Xf:  (6  n  ^C'K  61 ) 

1  »U15(E1  14)  ^PXN'<61)  »XMCCT{61) 

C-- . 

C -  input  data 

c - 

REA0C4i200)HMAXiITMAyilUPITEiI3UGi0ELMAX 

RPAD(4i2lO)XLiXLTiY?iCFL 

READ(4i220)FH01i='hC2iTl'T2iVOiVl 

READ(4i2'>2)!CHEMiT-,P0iyMA2 

REAn(4i223)P0iTHDiETAiTrLXiCF 

READ(4i?2l) ICONi IPT 

REAn(4,i?21)rCTrFFi:GFriINLETiIREADiITAPE 

RfAD(4i203)CAliyKCHEM 

CL0'E(UNlT  =  4ini??GG"  =  'FAVE*i£RR=99'5) 

I^dCTOFF  .EO.  P)  GO  TO  2C00 

REAP(78iP01)  (A(I)iU(Ii1)iU(Ii2)iL'(Ii3)i-C0V(I)iI  =  1iNMAX) 
CL05E(UNIT=78i0I6?0S'=’SAVE'iEPR=999) 


2000  CONTINUE 
C - 

C -  PRINT  INPUT  DATA 

C - 


MPIT6(fcli225) 

WPITE(61i?’r)MMAXiITMax,lwPITEiieUGiDFLMAX 

WRIT£(61i240)CFLaXLiXLTiYE 

WRITE(61i250)RHCliRhQ2iTliT2iVOiVl 

WRITE(61i252)ICHEyiTCiP0iyMA2 

WRITE{61i253)R0iTh0iETAiTPLXiCF 

WRITE (61  i2Fl  )  ICON  1  IFT 

WRITE(61i2f4)ICTDFFiir;E0ilNL£TiIREADiITAPE 
WPITE(61i256)CAliXKCHFM 
C - 

C -  NOZZLE  GEOMETRY 

C - 

A1  =  XLT 

IFdGEO  •E'-'.  0)  A1  =  XLT  -  XL 
DX  =  Al/FLnAT(NMAX-l) 

A2  =  XLT 

IFdGEO  .EO.  0)  A2  =  0. 

DO  5  I  =  liNMAX 

Xd)  =  XL  -  A2  ♦  OX*FLOATd-l) 

Ad)  =  <Y3  -  1.0)*X(I)**2/(XL**2)  +  1.0 
E  CONTINUE 

r. 

c -  INITIAL  CONDITIONP 

C - 

RU  =  S.BIE 
lEOPAT  s  0 


o  o  o 


111 


IV  V02.5  ^'on  :4-Jan-«3  13:46:37 


PAGE  00 


TIMF  =  0.0 
OT  *  0.0 
5W  =  O.n 
ITEP  =  0 
ISAVE  =  1 
cflmim  =  l.^=♦o^ 

NMIN  =  1 

IFdNLET  .>^0.  1)  N'‘iri  =  2 
00  10  I  =  l^N'IAX 

Udd)  =  fJHOl  -  PHD:*ry(l)  -  x(l))/  A1 
Ud»3)  =  T1  -  T2*(X(I)  -  X(l))/  A1 
Ud»2)  =  VO  V1*(X(I)  -  X(l))/  A1 
XMd)  =  Ud  .2)/':ORT(U(  d3)  > 

Ud»4)  =  0.0 
CNd)  =  0.0 
RXNd)  =  0.0 

10  CONTINUE 

IFdPEAO  .EO.  0)  GO  T"*  2010 

REA0(98»‘>C1)  (A(I).i;dd)»Ud<>2)»Ud»3)»Ud^4)»I  =  l^NMAX) 

CLnSE(UNIT=98»DISP0Sr=*SAVE'»EPR=999) 

2010  CONTINUE 

XMl  =  U(NM«X»l)*U(Nf  Ay»?)*A(N‘1AX) 

on  11  K  =  1  -.4 

FLUXS(K)  =  0. 

DO  11  I  =  l.NMAX 
Uld»K)=Ud»K) 

UlSd»K)  =  Ud»K) 

XNDOTd)  =  Ud  d  )  *Ud  *A  (I  ) /XNl 

11  CONTINUE 
RA2  =  PU/XMA2 
RO  =  RA2 

CVF  =  R0/(GAM0-1.) 


STAGNATION  CHENIS'^f’Y 


IFdCHEM  .FO.  1)  GO  TO  14 
PO  =  P0*1 .OF  +  05/ .9«t9 
RA  =  2.*RU/XNA2 
RA2  =  RU/XMA2 
CVF  =  3.*RA2 

Al  =  1.0E^06*RA2*T0*'’n«'GXP  {-THD/r0)/P0 
CAO  =  S0RT(A1  +  A1**:)/(A1  ♦  1.) 

RO  =  (1.  ♦  CAO)*RA^ 

RHOO  =  PO/ (PO*Tn* 1 .rE+06) 

GAMO  =  (RO  ♦  CVFI/CVF 
RDP  =  PO/RHOO 
Al  =  1.  -  CAO  ■ 

DSTG  =  (A1**A1)*(CAC**(2.*CA0)) 

IFdREAD  .ro.  1)  GC  TO  3 
00  13  I  =  l^NMAX 
A2  =  THn/(TO*Ud  »3)) 

IF(A2  .GT.  PS.)  Ud.4)  r  0.0 
IF(A2  .GT.  P5.)  GO  Tn  17 

Al  =  PO*EXF  (-THC/ (TC*IJ(  M3) ) )  /  (RH0C*Ud  il  )  > 
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f 

k 

< 

i 


IV  V02»5  ^’on  ?*-Jfn-S3  13546J57 


^AG?  003 


U(I»4)  *  (-A1  ■*■  SC9T(4,*A1  ♦  Al**2))/(2») 
IPdCTOFF  .cQ.  1)  U(!-*4)  =  CAl 

12  U1  { I  »4)  =  U(  N4) 

UIS(I»4)  =  U(r»i) 

9  =  RA2*(  1.  +  U1  <  I  i4  )  ) 
gam  s  (P  +  CVF)/CVF 

XM(I)  =  u  (  I  »2  ) /rOt?T  ( ♦UC I  »  3)  /  (GAMO^RO  > ) 

13  CCNTINUF 

IFdCHEH  .FC.  3  .0^.  Cd  .EQ,  0.0)  GO  TO  3 
U<1»4)  =  CAl 
Uld»4>  s  CAl 
UlSd»4)  =  CAl 


3  continue 

I - PRINT  INITIAL/PrUNDAPY  CONDITIONS 

C - 

14  WPITE(6l »r5F) 


WR1TE(61  ♦^bODX 

IFdCHFM  .ME.  1)WRITE(61»261)RA2»CVF»RH00»R0 
WRITE(61'»270)  d»Xd)»A(I)d=l»NMAX) 

WRITE  (61  ♦400)  ITEP»TIMF»r'T»IE09AT 
WRITE  (61  ♦410 

WR1TE(61^420  (I^X(I)^U(Ml)^Ud^2)VJ(I^3)^XMd)^ 

1  Ud^4)»Cf/(r)i0XN(I)^XM0CT(I)lI  =  l^NMAX) 

IP(I3UG  .:0.  0)  GO  TO  6C3 
WRITE(61^52C) 

WRITE(61^'ilO)(Uldd)»Uld^2)^Ul(I^3)^Ul(I^4)^I  =  l^NMAX) 
520  F0RMAT(//^1X ♦ 'MAIN  INITIAL  CONDITIONS  Ul’) 

60C  CONTINUE 

IFdCON.  EO.  0)  GO  TO  SCO 
WRIT£(61 ♦4’C) IPT 


300  CONTINUE 


ICOUNT  =  IWRITE 
ITEF  =  1 
100  CONTINUE 


C - 

c - FINITE  RATE  CHENISTPY 

C - 

TMIN  =  l.OE+06 
TMIN  =  XKCHEM*TNIN 


C - 

C - CALCULATE  TIME  fTFP 

C - 

CR  VMAX  =  0.0 

DO  60  T  =  1^NMAX 
P  s  (1.  ♦  U1(I^4))*PA: 

A1  »  (R**?  CVF*R)/ (CAfiO*RO*CVF) 

TEST  =  A9S(Ud^r))  SORT  (  A1*U(  I  ♦  3  ) ) 
IF(TEST  .GT.  VMAX)  VMAX  =  TEST 


o  o  o  o  o  r  *) 
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VN  IV  VO?. 5  Hon  ?4»Jsn->»3  13:46:37  >»AGE  004 

60  CONTINUE 

OT  =  CFL*nx/Vi1AX 
IFdCHFM  2)  CC  TO  62 

TFLCW  =  OX/VMflX 

62  00  61  I  =  l.NMAX 

R  =  (1.  ♦  U1 ( I »4 ) > ♦P A? 

A1  r  (R**2  +  CVF.R)/ (CA*'0*R0*CVF) 

CN(I)  =  DT*(AB£  (IJ  (  I  O)  )  +  £CRT(A1.U(  I  »3)  )  )/0X 
IC(CN(I)  .LT,  CFLMIN)  CFLHIN  =  CM(I) 

IF(CN(I)  ,CT.  1.0  ,r.r  .  CFL  .LT,  1.0)  GO  TO  63 
WRITE(61»440)  ITE^.I.OMCI) 

63  CONTINUE 

61  CONTINUE 

TIME  =  TIME  ♦  DT 


PREDICTOR  STEP 


ICALL  =  1 

00  70  I  =  liNMAX 

A?  = . THD/(T0*U1 ( 1.3)) 

IF(A2  .GT.  S5.)  CC(I)  =  0.0 
IF(A2  .GT.  85.)  GO  TO  70 

CF(I)  =  RO.tX*’ (-THn/(TC*Ul  ( I  .3)  ) )/ (  RH00*tJl  ( I  il )  ) 
70  CONTINUE 

IFdCHEM  .EC.  2)  CALL  THEM 

CALL  INTEG 

00  15  K  =  1.4 

00  15  1  =  NMIN  »NM AX 

Uld  .K)  =  Ud  .K)  ♦  DUd  .K) 

15  CONTINUE 


BOUNDAPY  CONDITIONS 


IPdNLET  .EG.  1)  G"  Tn  2020 

A1  =  RD*EXPC-THD/(T0*U1 (1 .3) ) ) /(RH00*U1 (1  .1 )  ) 

IFdCHEM  .ME.  2)  Ul(1.4)  =  (-A1  SQRT(4,*M  +  Al**2))/2. 
GAM  =  (4.  Uld  .4)  )  /3. 

COEF  *  GAMO/GAM 
A1  =  CAO  -  Uld. 4) 

A2  =  1.  -  Ul(1.4) 

A3  =  Uld. 4) 

A4  =  (RDR**A1)*(A3**(?,*A3))*(A2**A2)*EXP(A1) 

COEF  =  DSTn/A4 

Uld.l)  =  (COEF^Ul  (1  .3)**3)**(l./d.*A3)  ) 

Uld.2)  »  SQRT<2,*PA2*(4.-*-CA0-(4.  +  A3)*Uld.3) 

1  A1*THD/TO)/(3AMO*RO) ) 

2020  CONTINUE 

IFCIBUG  .EC,  0)  GO  T'j  610 
WRlTE(61.5no) 

WRITE(61.5lO)(Uld.l).L'ld.2).Uld.3).Uld.4).I  =  l.NKAX) 
500  FCRMAT(//»1X.»MAIN  f'f’EOICTGR  Ul') 

510  F0RMATdX.4E14.6) 

610  CONTINUE 


r>  o  o  o  o  o 


N  IV 


VO?. *5 


llA 

Mon  ?4-j3n-P3  13:46:37 


PAGE  COS 


C - CORPECTOR  STEP 

C-- - 

ICALL  =  2 
SW  =  1.  -  f.v; 

DO  76  I  = 

A2  =  THD/(Tn*Ul( I »  3)  ) 

IF{A2  .GT.  £'•-.)  CE(I)  =  0.0 
IF(A2  .r.T.  35.)  CO  TO  7f. 

CF(I)  =  Rn*FX  P  (-TH^  /  (  TC-.U1  ( I  »  3  ) )  )  /  (PHOO*U1  ( I  .  I  ) ) 
76  CONTINUE 

IFdCHEM  .EC.  2)  CALL  CHEM 
CALL  INTEG 
DO  24  I  =  l.NMAX 
RXN(I)  =  TCHEM(I) 

74  continue 

DO  25  K  =  1^4 
DO  25  I  =  NMIn.NMAX 

U1(I»K)  =  0.5*(U(I»K)  ♦  U1(1»K)  ♦  OU(I»K)) 

25  CONTINUE 


BOUNOArv  CONDITIONS 


IPdNLET  .EC.  1)  G'l  Tn  70?0 

A1  =  RD*rX“( -THP/ (  TO*i)  1  (  1  .3  )  )  )  /  {5H00*'J1  (1  .  1  )  ) 

IFdCHEM  .ME.  2)  Ul(1.4)  =  (-A1  ♦  S0RT{4,.A1  ♦  Al**2))/2. 
GAM=(4.'»-Uld»4))/3, 

COEF  =  CAMO/GAM 
A1  =  CAO  -  111(1.4) 

A2  =  1.  -  Ul(l»4) 

A3  =  U 1 ( 1 . 4 ) 

A4  =  {RC0**A1 ) *( A3** ( 2.*A3) ) • ( A2**A2 ) •£ XP (A  1 ) 

COEF  =  DSTG/A4 

Uld.l)  =  (CnEF*Ul  { 1  .  3  )  **3  )  *•  d  ./ (1  .fA?  )  ) 

Uld»2)  =  rCRT(r.*rAr*(4.+CA0-(4.^A3)*Uld  .3)  ♦ 

1  Al*THD/TO) / (GAM0*R0) ) 

203C  CONTINUE 

IFdBUG  .EO.  0)C0  TO  *20 
WRITE(61.530) 

530  FCRMAT(//.lXdMATN  CORRECTOR  Ul*) 

WRITE  (61. 510)  (Ul(I.l).l'l(I.2).Ul(I.3).Uld.4).I  =  l.NMAX) 

620  CONTINUE 


PRINT  OUTPUT 


IFdCON  .EO.  0)  GO  TO  3i 

WRITE(61.435)  ITER.U1(IPT»1).U1(IPT.2).U1(IPT.3).U1(IPT.4) 
GO  TO  35 

31  IFdTER  .NE.  ICPUNT)  CO  tq  35 
ICCUNT  =  lOO'JNT  ♦  IWPITE 
XMl  s  Ul  (NMAX  .1 )  *0  1  (*:*1AX  .2)  *A  (NMAX  ) 

DO  30  I  »  l.NMAX 
IFdCHEM  .e-Q.  1)  GO  TD  33 
P  =  PAr*( 1.  ♦  Ul ( I  .4  )  ) 

GAM  «  (P  ♦  CV^)/CyF 


C'>  o  o  o  o  o 
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\N  IV  V02.?  Hon  24"Jan-?3  13:46:37  PAGE  0C6 

XM<I)  =  U1  (  !  »?)/S0''T  (  r  ( I  ,  3) /(G«“0*P'3)  ) 

IFdCHEH  .FO,  2)  CT  ^7  3  2 
A1  =  U1 ( I »4) 

A2  =  TO^UK  J  n)/THr 

A3  =  («.+Al)*(l.+Al)*CAl*{U-Al)+3,*(2,-Al)eA2**2) 

A4  =  3.*A1*(1.-Al**r)*(l.  +  ?.*A.?) 

1  3.*(S.-*-3.*Al-Al**3)*A?**2 

A5  =  S0FT(A’/A4) 

XH(I)  =  A5*XM(n 
GO  TO  32 
33  CONTINUE 

XH(I)  =  U1  ( I »?)/SCrT (U1 ( I ,3) ) 

32  continue 

XMOOTd)  =  Uld»l)*Uld»2)*Ad)/XMl 
30  CONTINUE 

WRITE(61»400)  ITER»TIHE»DTiIEORAT 
WRITE(61»41C) 

WRITE<61»42C)  d»Xd)»Ul(I»l)»Ul(Ii2)»Ul(Id)»XM(I), 

1  Uld»4)»CNd)»'?XNd)»XMn0Td)d  =  l<NMAX) 

35  CONTINUE 


CONVEPPFNCE  TFf.T 


IFdTEP  .NE,  10)  GO  TO  1 5  <> 

SUM  S  O.c 
00  lAP  K  =  1»4 
DO  149  I  =  1»NHAX 

SUM  =  SUM  +(<U1(I»K)  -  UlSd»K))**?) 
149  CONTINUE 

SUMl  =  SORT(SUM/FLnAT(N'‘'AX)) 
lf9  CONTINUE 

IFdSAVE  .NE.  10)  Gr  Tn  41 

ISAVE  =  0 

DEL  s  0.0 

SUM  s  0.0 

DO  40  K  =  1.4 

00  40  I  =  l.NMAX 

TEST  =  ABS(U1(T.K)  -  UlSd.K)) 

IF(TEST  .GT.  DEL)  OFL  =  TEST 

SUM  =  SUM  ♦((Uld.K)  -  llSd.K))**?) 

40  CONTINUE 

SUM  s  SOPT(SUM/FLO«T(NMf X)) 

SUM  s  AL0G1C(SUH/SUM1) 

IFdTAPE  .PC,  2)  WPrTF(77,900)ITEP»SUH 
IF(DEL  .LT.  PFLHAX)  F.n  TO  99S 
00  42  K  =  1»4 
DO  42  I  =  l.NMAV 
UlSd  .K)  =  Ul  (I  .K  ) 

4?  CONTINUE 

41  CONTINUE 


RESE'^  AFtlAYS  ‘^n  nFu  VALUES 


ITEP  s  ITFP  ♦  1 


ooo  ooooo 
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ISAVE  =  IfwIVE  •»  1 

IFdTEP  .GT.  UMAX)  3-  999 

00  45  K  =  1»4 

00  45  I  =  liNMAX 

U(I»K)  =  UKIiK) 

45  CONTINUF 
GO  TO  lOP 

r 

END  PF  LOOP  C 

C 


9‘'8  XMl  =  U1  (NMAX  d  )  *11 1  (►;►' AX  »?)  ♦  A  (NMAX  > 

DO  55  I  =  1»NMAX 
IFdCHFM  .EO,  1)  GP  ^  3 
R  =  RA?*d.  *■  U1  (  I  -(A  )  ) 

gam  =  (f  ♦  cyF)/cvp 

XMd)  =  UK  I  12)/S0PT  (';flM*R*ul  d  *  3)  /  (GAM0*R0)  ) 

IF(TCHFM  .'^0,  2)  GO  T'*  54 
A  1  =  U 1  (  I  1  4  ) 

A2  =  TO^Ul ( r ♦3)/THC 

A3  =  C4.*'Al)*(l,t-Al)*(ri*(l,-Al)-*-3.*(2.-Al)*A?**2) 

A4  =  3.*A1*(1,-A1**:)*(1.+?.*A2) 

1  +  3,*(P,-<-3.*d-Al**3)*A2**? 

AS  =  SCRT(A3/A4) 

XHd  )  =  A5*XM(  I) 

GO  TO  54 
5?  CONTINUE 

XMd)  =  U1  (I  »2)/5n9T  (UK  I  »3)  ) 

54  continue 

XMCOTd)  =  UKI»l)*UKI->2)*Ad)/XMl 

55  CONTINUE 

WRITE(61»4C0)  KEFdlMF^DTdElRAT 
WRITE(61 »41C) 

WRITE(61»4?0>  d»Xd)»UKI»l)»UKI»2)»UKI»3)»XMd)» 

1  UKI-*4)iCN(I)»RXNd)»XMD0T{I)»I  =  KNMAX) 

999  CONTINUE 

WRITEC6K445)  CFLMlf: 

CL0£ECUNIT  =  61»DICP!:r.E  =  '5AVE') 

IFdTAPE  .PO.  l)WRITr(9O>901)(Ad)^UKId)»UKI*2)*UKK3)» 
1  UKI  *4)  d  =  KNMAX) 

STOP 


FORMAT  statements 


2CC  FCRMAT<4I5  »F 10  .G) 

205  FnRMAT(2=’10.0) 

210  F0RMAT<4F10.C) 

2?0  F0RMAT(6F10.0) 

nn  F0RMATC5I5) 

221  F0RMAT(?I5) 

222  FrRMATd5»3FlO.O) 

223  FORMATOr  lO.OdElP.I) 

224  FCRMAT{4E14.6> 


IV 


V02.5 


Mon  24-J--;n-‘<3  13:46:37 


PAGE  008 


225  PCPfMT(//»lX»*  f'ATA**//) 

230  FPR^lATdy^'NMAXs'^I-OX^'IT^AXs'^IcOX^'IWRITEa'^ISOXt 
I  •I?UG='»I5i3Xi'D5L»-'AX=*»810.3) 

240  FrR»AT(/»lX»*C=’L='iFlC,3»3X^*XL=*dl0.3»3y, 

I  •XLT=»»Eir',-'»?X»*yc  =  *  »E10.3> 

250  FOr?r;AT(/»lX»'FHCl='i'’18.10»3X»'RH02=dElC.3»3X»*Tl  =  »»E18.10»3 

1  •T2=*  »E10,3i3X»‘V0  =  *  »El8.10»3XdVl=' 1810.3) 

251  F0RMAT(/ilXi*IC0N='iI5i2Xi'IPT=*»I5) 

254  F0RMAT(/ilXi'IC'^CFF  =  'irEi3Xi*IGE0='ir5i3Xi*INLET  =  ‘»I5i3X» 

1  • I  READS’ 1 1  5 1 ?X  1  • ITAPE=  •  il 5) 

252  F0RPAT(/i1Xi’ICH'M='iI5>3Xi’T0='iE10.3i3Xi'P0='iE10.3» 

1  3X  »’XMA2=  •  I'^IO.B) 

253  F0RMAT(/ilXi*PD=’iriC.3»3Xi’THDs'i£10.3i3Xi»ETA=’iE10.3» 

1  3X»’TPLX='iFl0,3»3X»’CF=’»E10.3) 

255  FPRMAT(//»’ INITIAL/“CUNDAPY  CONO 1 T 1 ONS ’  i / / ) 

256  FCRMAT(/ilX»'CAl  =  ’ iElC.3i3Xi 

1  'XKCHEMs’  »E10.3) 

260  FPRMATdXi  *DX=’ iFin,3i/) 

261  FORMAT  (/«1X*’PA2='  ■.']®.lO»3Xi'CVF=’i£13.lOi3X» 

I  ’RHaOs’iFn.lOiBXi’ROs'iEl^.lO) 

270  FGRMATdXi’I  =  'iI3i?Xi’X='iE14,6i3Xi'A  =  ’iE14.6) 

4P0  FCRMAT(//»lXi«  IT"'^=’ iI4»3Xi ’TIMEa’ ♦El4.6i3X»’0T=' icl4,6 
1  i3Xi* lEQPATa'  iI4) 

410  F0RMAT(/i4Xi'T’i8Xi'X’il2Xi'PH3'd2X»’U’»13Xi'T'il3X»'M'»12X» 
1  •CA'il2X»’CN’d2X-<'RXN’il0Xi’XMi)0T') 

42C  FPRMAT(3XiI3iRE14,6> 

430  FCRMAT(//ilX»'CCr.’VEPGENCF  C  NECK -U 1  »  U2  ♦  U  3  •  »  3X  »  •  I  a  ’  »  I  5  »  /  /  ) 

435  FpRMATdXiI5i3Xi4T14.6) 

440  FGRMAT(//ilXi'ITEr='iT4»3Xi’I='iI4i3X»'CM=’»E14.6i//) 

44E  F0RMAT(//»1Xi*CFLM!N=',e14.6) 

441  FCRMAT(//ilXi’IT-F=*iI4i3Xi’I=’iI4i3Xi’RXNa’,El4.6»//) 

ROC  FCRNATdXiI5»ri4.6) 

ROl  F0RFAT(5F14.6) 
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SUPROUTINF  INTPC 

COMMON  Ul(E!»4)»«)((4»4),lXIf/V(4<4)»P-j(6l^4)»U(61»4)»A(ol)» 

1  C»rFL»rT^rXfnM'0»II»N’^AXiSW»I5UG»PLUXS(4)MCALL 

2  iT-jnD 

3  »0MECA(61)iTCHrM(61)»T(61>»ACHEM»CVF»F.TA»ICHEM 

4  ♦'=A2»PD^PH0n»^C^TH')»TRLX»TC»CE(61)»nCD91»DCDTlt0MP»I 

5  »ECr’N(61)»ICTnPr»XKCHIM»INLFT»NMIN 
DIMENSION  DA(4)»0AINV(4)»FLUX(4)»W(4)»X(4)»Y(4)»2(4) 

C - EXPLICIT  STAGE 

C - 

IFdCTOFF  ,EQ,  0)  GC  TO  2000 
lEOPAT  =  0 


DO  80  I  =  2»NMAX 

IF(I  ,NE.  NMAX)  GO  TO  71 

C4  =  U1 ( I »4)  -  U1 { I-l »4) 

GO  TO  72 

71  C4  =  (1.-SW)*U1<I  +  144)-(1,-2,*SM)*‘J1(I»4)-SW*U1<1-1»4) 

72  0U<I^4)  =  -nT*(Ul<  I  ^2)  *C4) /DX  OT*OMEGA(I) 

OEQ  =  EfONI I )  -  U1  (  I  '4) 

IP(AnS(r!U(I»4))  ,GE.  ArF(DFQ))  lEORAT  =  I 

80  continue 
2000  CONTINIJ- 
ISND  =  0 
DCPPl  =  0.0 
PCOTl  =  0.0 
00  ?  I  =NMIN»NMAX 
IF(I  .NE.  1)  GO  TO  ? 

Cl  =  U1<I+1»1)  “  U1(I»1) 

C2  =  U1(I  +  1»2)  -  UKI-*:) 

C3  =  Ul(I•^l13)  -  UK  I-.2) 

C4  =  UKKK4)  -  Ul(!»4) 

CP  =  ALOG(A( I+l) )  -  ALrC(A(I)) 

GO  TO  4 

2  IF(I  .ME.  NMAX)  GO  TO  3 
Cl  =  UK  IK)  -  UK  I-l  K) 

C2  =  UKI»?)  -  UKI-K2) 

C3  =  UKI»3)  -  UKI-K?) 

C4  =  UKI»4>  -  UKI-1»4) 

C5  =  ALOG(A(I))  -  ALnC(A(T-l)) 

GO  TO  4 

3  continue 

Cl  =  (l.-SW)*UKKlil)  -  <1.-2.*SW)*UKI  »1)  -  SW^UKI-IK) 

C2  =  (l.-SW)*UK  K1  »2)  -  (l.-2.*SM)*UKK2)  -  SW*UKI-1»2) 

C3  =  (l.-SW)*UKI  +  K3)  -  (l.-2.*SW)*UKTi3)  -  SW*UKI-1»3) 

C4  *  (l.-SW)*UKKl»4)  -  (l.-2.*SW)*UKI  »4)  -  SW*UKI-1»4) 

C5  =  <1.-SW)*ALPG(A(K1))  -  (  1 .  -  2  .  •  SN  >  •  ALOC.  (  A  (  I )  ) 

1  -  SW*AL06<A(I-1) ) 

4  CONTINUE 
OCOP  =  0.0 
DCDT  8  0.0 
DU(I»4)  =  0.0 

IFCICHFM  .EO.  3)  r-n  T  ■>  04 
IFdCHFM  .EC.  1)  OMEOAd)  a  0.0 


o  o  o 
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0U(I»4)  =  -0T*(U1(I»2)»'C4)/DX  ♦  OT-OMFGAd) 

ICUT  =  0 

IFdCTOFF  .£n.  0)  C"  ‘"H  55 

OFO  =  XKCHrM*(F.CON(I  )  -  J  1  (  I  i  4  )  ) 

CON  a  U1 ( 1 14 ) ♦*?/( 1 .-Ul (  I  »4)  ) 

CONP  a  U1 ( T  »4)  *  DU( r  »4) 

EQCPIT  =  U1(I»4)  +  TFn 

IF(CONP  .LP.  O.C  .CP.  CCNP  .GE.  1.0)  G3  TO  58 
EQPAT  =  CQN/CE(I) 

IF(I  .GT.  TFQPAT)  GO  TO  55 
58  ICUT  a  1 

0U(I»4)  =  OEQ 

QMEGA(I)  =  DU(I»4)/CT  +  U 1 ( I » 2 ) aCA/DX 
Q(I)  =  OMFGA(  DaTHP/ (  ^.♦TO) 

55  IFdCALL  .EQ.  2)  GO  TC  56 
TCHEPd)  =  1.0 

.  IFdCUT  .EC.  1)  TCHFMd)  =  3. 

GO  TO  60 

56  IFdCUT  .EO.  0)  GO  TO  60 
IF(TCHEHd)  ,E^.  3.0)  CO  To  57 
TCHEMd)  =  2.0 

GO  TO  60 

57  TCHEMd)  =  4.0 
GQ  TO  60 

C - 

C - EQUILI6FIUM 

C - 

5*  A1  a  CFd) 

A2  a  U1 d  44 ) 

Uld44)  a  (-A1  '■•0FT(4,*A1  -f  41*42))/2. 

nu(l»4)  a  Uld44)  -  S2 
OMECACI)  a  0.0 

Q(!)  a  0.0 

C4  a  0.0 
Al  a  U1 (I  44) 

A2  a  THP/TO 

A3  a  Al*d.-Al)/(?.-M) 

DCOP  a  -A3/Uld4l) 

OCOT  a  A2*A3/(Uld43)**-’) 

IPCI  .EO.  1)  PCCRl  a  nccf? 

IFd  .EO.  1)  DCOTl  a  OCGT 
TCHF.Md)  a  4.0 


NONEOUILI9RIU?’ 


60  continue 

R  a  RA2ad.  ♦  U1  d  44  )  ) 

31  a  l./(l.  ♦  TH0*''C0T/(  3.*T0) ) 

3?  a  THD*0C0R/<?.*TD) 

0U<l4l)  a  -OTaCin  {  !  42  )  4C1  ■»  Ul<l4l)*C2  4  U 1  d  4  1 )  *U  1  d  4  2  )  aCS  )  /  D  X 
CUd42)  a  -CTaUl  ( I  42)  aC’/DX  -  OT*  (  R  aC  laR  A2  *0  1  d  4  3  )  •  (  C4 

1  ♦OCCFaCl  a 

2  ♦  PaUl(T4?)aCl/Ul(I 4l))/(GAM0aR0a0X) 

IP{ICHEH  .FO.  1)  r(I)  =  0.0 


o  o  o  o  cT  r>  o  o 
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OU(IO)  =  -PT*U1  (  1  O)  •C3/rx  -  DT*T1*C7*U1(I  .3)/CVF  - 
1  •(C?  ■*•  U1  (I  i2)^C5)/DX  -  DT*2(I) 

5  C0NTIf'lJ=- 
0U(1»4)  =  ^.0 
IF(I5UG  .F'^.  0)  GO  TO  2^0 
WRrTE(61 »ior) 

100  FORMAT (// »ix »• iNTfG  EXPLICIT  DU*) 

WRITE(fel»110)(DU(I.l).DU(I»2)iDU(I»3)OU(I»4)iOMECA(I).CE(I)» 
1  I=liNMAX> 

110  format ( IX .6E14  ,  D  ) 

200  CONTINUE 

IFdNLFT  ,^r.,  C)  CALL  tmory 


STARTLITY  CHECK 


IF(CFL  .LF,  1,0  PFTUPr 


BOUNOARY  CONDITIONS 


nn  47  K  =  1.4 

FLUX(K)  =  FLUXS(K) 

47  CONTINUE 
ISAVE  =  2 

IFCINLET  .ro.  O  CD  '"0  2010 
ISAVE  =  0 
DO  410  K  =  1  »4 
FLUXS(K)  =  0.0 
41R  CONTINUE 
*'010  CONTINUE 

IF<£W  .FQ.  1.0  .ANf'.  ■'C"LL  .  E  ) .  1)  ISAVE  =  NMAX-l 
IF(FW  ,FR.  O.C  .AND.  TCALL  ,E0,  ?)  ISAVE  =  NHAX-l 


IMPLICIT  STAGE 


IBND  =  1 

DO  40  N  rNMlN.NMAX 

I  =  SU*N  (l.-rW)*(‘IMAy-t>NMIN-N) 

II  =  I 

IF(I  .NE,  1)  GO  TO  13 
CALL  PNDRY 
DO  11  K  =  1.4 
X(K)  =  0.0 
00  11  L  =  1.4 

X(K)  =  y(K)  ♦  Sy(K.L)«FLUX(L) 

11  CONTINUE 

DO  12  K  =  1.4 
FLUX(K)  =  X(K) 

1?  CONTINUE 
13  CONTINUF 
CALL  EVECT 
00  10  K  =  1.4 

W(K)  =  OU(l.K)  f  DT*FLIJV  (K)/DX 
1C  rCNTiMU'" 

00  IS  K  »  1»4 
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X(K)  =  0.0 
00  15  L  =  1»4 

X(K)  =  X(K)  SX(K»L)*W(L) 

15  CONTIMJE 

0A(1)  =  APS(Ul(I»r)  -  fCHEM*C)  -  (DX/OT) 

0A(2>  *  AariUKI'?))  -  (OX/CT) 

0.^(4)  =  AP*' (U1  (I  ^?)  5  -  (OX/DT) 

0A<3)  =  ACSCUKIi?)  ♦  ArHEH*C)  -  (DX/OT) 

IPdCHt^M  .NT.  2  .Of’.  T  .EO.  1)  nA(4.)  =  0.0 

DO  ?0  K  =  1»4 

OAIK)  =  AMAX1(0A(K) »0.C) 

DAIKV(K)  =  l./(l.  ♦  nT*DA(K)/OX) 

IPCICHE'^  .K'E.  2  .OR.  I  ,E0.  1)  0AINV(4)  =  0.0 
Y(K)  =  DAInv(K>*X(K) 

20  CCNTINUF 

IFCI  .NE.  !)  GC  TO  24 
A1  =  DAINV(I) 

A2  *  ((U1(I»2)  -  ACHE“*C)  /  (U1  ( 1 12)  ACHEM  *C  )  )  *0  A  (  1 )  •  A 1  *01 /O  X 

Y(l)  =  AI«X<1  ) 

Y(?)  =  X(2> 

Y(3)  =  Xd')  -  A2*X(1) 

Y(4)  =  C.O 

24  or  25  K  =  1.4 
0U<1.K)  =  0.0 
on  25  L  =  1.4 

OlJd.K)  =  'UH.K)  .■  5X1MV(K.L)*Y(L) 

25  CONTINU- 

00  30  K  =  1.4 

2<K)  =  OA(K)*Y(K) 

30  CONTI  HUE 
00  35  K  =  1.4 
FLUX(K)  =  f'.O 
00  35  L  =  1.4 

FLUX(K)  =  FLUX(K)  ♦  SX  INV  (K  .L  )  *2  (D 
35  continue 

IFd  .NE.  ISAVE)  Cr  TO  ?7 
00  38  K  =  1.4 
FLIIXS(K)  =  FLUXCK) 

38  CONTINUE 
37  CONTINUE 

IF (I  BUG  .E”.  0)  GO  TT  210 
WRITE  (61 .120  I  .f'X  .OT 

12C  FCRMAT(//.1X.' INT'^C  I-i“LICIT  STE  P  •  »  3X  .  •  I  =  *  .  1  3 . 3X  .  •  0X=  *  . 

1  F7,4.3X»'2T='.F7.4./.1X.»W.X.Y.Z.0U»FLUX') 

WRITE(61»130{W(K).X(E).Y(K).Z(K).DUd.K).FLUX(K)»K=1.4) 

130  FCRHAT( IX.6E14.6) 

210  CONTINUE 
40  CONTINUE 
RETURN 
FNP 
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SU"*:  nuTI‘J>  EVECT 

COMKON  U1(61»4)<EX(4»4)iEXINV(4<4)0!J(61»4)»U(^1»4)» 

1  4(61)»fiCFL?''T^TX?GA‘^.'!iII»h'i\X»3UiI2U^-»rLUXS(4)»IC9LL 

:  »  I  F  N  c 

3  ^OMTCiA  (£  1  )  »TCH^‘  (  6  1  )  O  <  6  1  )  1  ACHE*^  »C  Vr  »ETA  »  irH£*1» 

4  RA.?»F,  nirHC0^rr»THC»TRLX»T0»CE(6n‘i?C0K1^0C0Tl^C'lP»IEC 

F  »FC0N(61  )  ^IC'^rr=-,XKCHE'1 

I  =  II 

C  =  SORT (U1 ( I »3)  ) 

R  =  RA2*  (  1  .  ♦  U1  (  1  •*4  )  ) 

OCDP  =  C«0 
nCDT  =  0.0 

ACHEM  =  £0'’T(  (CVF*"  ♦  R  *  ♦  2  )  /  (  G  0  ♦  R  0  »  C  VF  )  ) 

IFdCHEM  ,MF,  3)  GP  TG  ? 

A1  =  U1  (  I  <4  ) 

A2  =  T0*U1  ( I  -*3)  /THD 

A3  =  GAMG^^C*  C  A1  *  (  1  .-Al  ) -^3.  •  (  ?,-Al )  *  A2**2  ) 

A4  =  A1*(1.-A1**?)*(  1  ,i-2.*A2)  +  (8.i-3.*Al-Al**3)*A2**2 
ACHFK  =  50FT(RA2*A4/A3) 

B1  =  U1  (  I  »4  )  •(  1  .-Ul  (  M4)  )  /  (2.-U1  (  I  »4)  ) 

C2  =  THP/TO 

OCPR  =  -R1/L'1(I»1) 

DCOT  =  C2*°1/(U1(I' 

•’  CONTINUE 

A1  =  l./(l.  +  THP^RCrT/O.^TO)) 

A2  =  THP*OrrR/(’.*TC  ) 

IFdCHEM  ,FC.  1)  =  1, 

IFd?L'G  .EC.  0)  r.P  TO  IP'I 
URdE  (61  »  ICC)  I  »  ACHFV  ■.C 

10  0  FrRMAT(//dX»'EVrfT'dX»d  =  'd3dX-*'A  =  *»E14.6»3XdC=*^El4.6) 
200  CONTINUE 

C -  CALCULATE  MATRIX  OF  EIGENVECTGPS 

C - 

SXd»<^)  =  c. 

SX(?»4)  = 

SX(3»4)  =  C. 

5X(4«1)  =  0. 

SX(4»2)  =  C. 

SX(4»3)  =  0. 

SX(4»4)  =  C. 

SX(1  »1 )  =  1  . 

SXd»2)  =  -ACHEM*C*G  AMC*R0*U1  (I  d  ) /(R*U1  (  I  »3  )  ♦ 

1  FA2*U1  ( !  d)  •Ul  (I  »1)  •POOR) 

SX(1»3)  =  (Uld»l)/Lld»3)) 

1  *(R4.RA?*Uld»?)*PC0T)/(R^PA2*Uld»l)*0CDR) 

SX(2»1)  *  1. 

SX<2»2)  =  C. 

SX(2»3)  =  -L’l  (I  ♦  1 ) /(  A  !  *(0*01  d  »3) /CVF  -  A2*Uldyl)}) 

SX(3yl)  =  1. 

SX(?y2)  =  -«X(ly?) 

SX(3i3)  =  •■XdO) 

SX(4yl)  = 

SX(4»2>  = 


•fc-r  >■ 


vo2.‘; 


ron  r4“jsn-83  13:49:13 


“AGE  OC? 


SX(4,3)  =  0, 

IrdCHEM  .Mr.  2  .09.  r  .=0.  1)  G?  TO  10 

SX  (1  »4)  =  ?*U1  (  I  »1  )  /p 

SX(?»4)  =  0. 

SX{3  »4)  =  <-X  (  1 ,4  ) 

SX<4»4)  =  1. 


c - 

C - CALCUL-'Tr  irvr'’S'  f’aTRlX 

ir  A1  =  fX('»2)*SX(ld)-5X(ld)»PX(’»3) 
1  -5X(2»?)*(SX(?.2)-SX(li?)) 


SXlNVd^'^)  =  0. 

SXInV(2»4)  =  n. 

SXIKV(?»4)  =  C. 

SXInV(4,1)  =  0. 

sxrNV(*»2>  =  0. 

SXINV(4»3)  =  0. 

SXlNV(4»A)  =  0. 

SX!^'V(l-.l)  =  -SX(3»2)*‘'y(2d)/Al 

SXiNV/di?)  =  (SX(3»2)*:XdO)-SX(  Idl^SXOdll/Al 

SXlMVdd)  =  SX(d2)vPx(2i3)/Al 
SXlMV(2d)  =  (SV(?d)  -  SX(3i3))/Al 

SXlNV(2d)  =  0. 

SXlHV<2d)  =  -SXrMVCOd) 

SXTNVOtl)  =  SX(’i2)/d 

SXINV(3^2)  =  (SXCld)  -  SX(3»?))/A1 
SXIh'V(?<3)  =  SXTNV(?»1) 

IFdCHE^'  .VE.  :  .CP.  T  .EO.  1)  GP.  TO  20 
A2  =  l./(P  ♦  CVD 

SX1NV(1»4)  =  -A?*rA2»CvP*Ul(I»l)/P 
SXINV(?»4)  =  0. 

£XINV(3i<*)  =  -Aro°A2*Ml  (I  i3) 

SXlNV(4d)  =  0. 

S  X  I N  V  ( 4 . 2  )  =  0  . 

SyiNV(4»3>  =  0. 

SXINV(4,4)  =  1. 

20  CONTTHUE 

IFdaUG  .EO.  0)  GC  TP  210 
4PITE(6]  dlO) 

IIP  F0PMAT(//>lXdEVECT  '  X  .  S  X  I  N  \/ •  ) 

WRITE(61d20)(SX(Ld)»:x(H2)»SX(H3)»SX(L^4)iSXlNV(Ld)» 
1  r,X  If'V  (L  .2)»?XI\V(L»3)<SXINV(Li4)»L=1»4) 

120  format  dX  t8E14.6> 

210  CP'NTINUE 
RETURN 


UJ 


12A 


•i  IV  V02.''  Vcn  :4-Jan-??  n:49:r3  PAGE  OOl 

SU9P0UTI\r  ■=N')PY 

common  Ul(6li4)»9X(4^4)»SXlNV(4»4)»DU(bl»4)«ll(61»4)i 
1  A(61)^C»C-LtCT»DX»GA'lJiIIiNM«X»5W»!2UG»FL',;XS(4)»ICALL 

t  ir  ND 

»^M^GA(tl)^TCHrM(tl)^^(61)^ACHEM^CVF^tTA^ICHEM^ 
P"?»PD»'’H0n-.'"C^THr.»T'’L':»T0^C2(61)»DCDPl»DCC'Tl»0MP^IEGr.' 
»rCGh(61)»ICTn>rr,XKCME'1 
OIMENSIDn  GU1{4) 


CALCULATE  P1P2  .-'ATPIX 


R  =  (1.  ♦  L'1(1»4))*FA2 

Z  =  Ul(1^4) 

A1  =  SCFT((P**2  ♦  CVF*n ) / ( GAM0*R0*CVF  )  ) 

IFdCHEM  .NP,  3)  GC  Tl  ? 

01  =  U1 ( 1  14) 

B?  =  T0*Ul(l»'»)/ThD 

B3  =  GAMO*bo*{P,i«(1.-"1)  +  3.*(2.-B1)*3’**2) 

B4  =  El*(l.-Bl**2)*(l.i-?.*P2)  + 

1  (8.  +  3.*r!-£l**3)<‘-3  2**2 

A1  =  SCPT(,7A?*B4/^3) 

2  CCNTIMUF 

A2  =  (UKliT)  +  TH'/T’"')  cOCPPl 
AX  s  (UKld)  +  THr/TOvDCPTl 
A4  =  GArO*PC/RA2 
AS  =  THri*DCCPl/(T:  *'0  1  ( 1  *3)  ) 

A6  =  THP*DCnT 1/(T0*U1 (  1 13)  ) 

A7  =  A1*GAHO*=!0«U1  (  1  1 1)  *SPPT(U1  (  1  1  3)  )  /  {R*U1(  1  i3) 

1  ♦  P  A2«'l!l  (  1  »  3  )  *IJ  !  (  1  1 1  )  *PC0R  1  ) 

A9={Ul(l»l)/Ul(li?))*d*RA2*Ul(l»3)*DC3Tl)/(Pi-RA2*Ul(l»l)*r)CCPl) 
A8  =  (-(  1  .-t-Z  ) /U1  ( 1  1 1  ) +  A5  )  ♦  (- A7*  (  4  ,i.Zi.A3  )  -  A  4  »  U  1  ( 1  i  2  )  *  A9  ) 

1  ♦  (3./L'1(1i3)i-A6)*(A4»J1(1,2) 

2  ♦A2*A7) 

Pll  =  A4*Ul(li2)*C>,/'Jl(l»3)+A6) 

P12  =  -A2*(3./U1  ( 1  i3)  i-f  6  )  +  (-(  1  .i-Z) /U1  ( 1  »n +A5)* 

1  (4.+2+A3) 

P17  =  -A4*Ul(l»2)*(-(l,  +  Z)/Ul(l»l)i-A5) 

SX(1 »1)  =  Fll/A? 

SX<1»2)  =  -P11*A7/AS 
SX(li3)  =  P11*A9/Afi 
SX(1»4)  =  0.0 
SX(?»1)  =  P12/A? 

SX(2»2)  =  -P12*A7/AR 
SX(2*3)  =  B12*A9/A? 
sx(?»4)  =  n.o 
SX(3»1)  =  P13/A« 

SX(3»2)  =  -P13*A7//'P 
SX<3»3)  =  P13*A9/A? 

SX(3t4)  *  0. 

SX(4»1)  =  C. 

SX(4i2)  =  C. 

SX<4,3)  = 

SX(4»4)  =  C. 

IFCIBNp  .EO.  DrETliFU 


I 


o  o  o 
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S  IV  VC2.5  Mon  24-J.'n-??  13:49:53 


PAGE  C02 


CALCULATE  Dl'(l»K) 


DO  10  K  =  1  i4 
DUKK)  =  0. 

DP  10  L  =  1 i4 

OUl(K)  =  nuKK)  SX  (K  )  *00  (  1  »L  ) 

10  CONTINUE 

00  20  K  =  li4 
D  U  (  1  » K  )  =  n  u  1  ( K  ) 

20  continue; 

IP(I9UG  .EO.  0)  GO  TO  ?G0 
y9ITE<61  ^If^P) 

Wt?rTE(61*110)(£X(K:,l)'rx(K»?)»5X(K»?)^SX<K»4)»0U(l»K)»K  =  l»4) 
UPITE(61»12C)  Al»Ar»A3»A4,a5»A6iA7»A3»A9 
100  FOP“AT  ( // »1X  »  •  3r:0RY  P1P2»DU1') 

Ilf  FfPMATC  l>:.5E14.h) 

12C  FCRMAK/slXiOrit.c) 

20C  CONTINUE 

oeturn 

END 


o  O 
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'I  IV  VO?.*;  Mon  ?4-jf.n-n3  n:so:30 


PAGE  001 


SUPPOUTIN'  CH'^M 

common  lll(frl»4)»«X(4»4)«SXTMV(4i4)»PU<6  1»4)«U(61*4)» 

1  A  (  6  1  )  »  r  »  :  -  L  1 »  D  X  <  G  A“  0  ,  I  I  ,  rr-l  A  X  1  -i »  I  :u  o  »  F  L  U  X  s  (  4  )  1  I  C  A  L  L 

2  I'^,01^^•^GA(61)»TCHE“(61)^0(61)^ACHE(M^VF»cTA^ICHf.  Ml 

3  RA?»Pr^RHOr»PO^THf'»TOLX»TO»C=(61)OCD‘?l»OCDTlOMp^IEC 

S  »FCCN(61)»!CTCFF^XKCHEri 


CALCULATE  ChEMICAL  PRODUCTION  RATE 


00  10  I  =  l^MMAX 
IF(CE(I)  .►'F.  O.C)  0  0  TO  2 
OMECAd)  =  C.O 
13(1)  =  0.0 
GO  TO  1C 

2  A1  =  EX?(-THD/(T?<:U1(!»3))) 

A2  =  UlCIi4)*«2/<l.-lil(T  ♦4)) 

A3  =  A2/CZ ( ! ) 

CMEGA(I)  =  ( 1  .-A3)  ♦  (1  .-I'l  (  T  <4)  )*A1  ♦TRLX*U1  (  I  il  )*U1  ( I  »3)**ETA 
0(1)  =  OMFOAd)  ♦TH  [1/(3. ♦TO) 

10  CONTINUE 

IFdPUG  .F-.  C)  GO  TO  i’f'O 
WRITE  (61  dOP) 

UP  ITF  (61  dl  0  )  (L’l  (I  d  )  ^L'l  d  »3  )  I'Jl  (  I  *4)  O.-IFGA  (I  )  d  =  l  »NMAX) 

20C  FORMAT  (// d  CHF'^  A  i  T  C  HF  M  •  ) 

210  FnPHATdXdFlP.lC) 

220  FCRMAT(iyd4dtl4.6) 

3C0  CONTINUE 
’FTl'RN 
END 


fll 


